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METHOD OF MULT I -DIMENSIONAL ANALYSIS OF VISCOELASTIC 
MATERIALS FOR STRESS, STRAIN, AND DEFORMATION 

FIELD OF THE INVENTION 

This invention relates to the methods and processes of 
5 multi-dimensional analysis of material properties, and the 

methods and processes of solving any arbitrary multi-dimensional 
scientific or engineering design problem requiring solutions for 
stress, strain and deformation which demand incorporation of a 
material 'constitutive function' into the solution, wherein the 

10 solution for that constitutive function, which defines the 

p 

\Q relationship between stress and strain in solids, and stress and 

|* rate of deformation in fluids, incorporates tensor valued 

in 

y coefficients and scalar valued constitutive functions, where the 

ly , ... 

ig tensor valued coefficients and the constitutive functions are 

15 I*j material specific. 

m 

*D BACKGROUND OF THE INVENTION 

§~i 

^ The Need for the Higher Order Multi-dimensional Constitutive 
Functions of Viscoelastic Materials 

Hooke's law of linear elasticity and the Navier-Poisson 
20 equations of Newtonian fluids represent idealized limits of the 
more generalized and more complex theory of material behavior 
known as ' viscoelasticity ' . 

Viscoelastic constitutive functions are time dependent, 
where the one dimensional relationships between static stress a 

1 

\ 



and a resulting strain rate €(t), or a static strain e and the 
internal state of stress a(t) are, 

o(t) - f(t,e) e(t) - £*(t,o) 

Idealized one dimensional linear models for these materials were 
first introduced in the last half of the 19th century by Maxwell, 
Meyer, Kelvin, Voigt, and others (c. 1868 - 1889) . 

The most generalized form of a three dimensional linear 
constitutive equation for static loading or deformation may be 
written as, 

°<*>ii " fijkiltki' e w ) , elt)^. - f* ijkl (t kl , o kl ) 

where the f * <t kl ,a kl ) ijkl and f (t kl , e kl ) ijkl are highly non-linear 
functions, generally assumed to be independent. However, these 
functions can not be independent, because of chemical bonding and 
the number of independent functions associated with various 
degrees of material symmetry have not been determined. Because 
multi-dimensional behavior is not well understood theoretically, 
scientific and engineering applications for these types of 
materials have been constrained and the characterization of 
material properties more of an art than a science. 

Experimental work conducted on a limited number of linear 
viscoelastic solids has shown that the matrix of functions, the 
f*(t kl ). jkl and f(t kl ) ijkl , may occasionally be approximated by a 



matrix of independent coefficients and a single scalar valued 
creep compliance or relaxation modulus function, 

f* M m -> f*(t) , f(t) ijkl -> td^lfU) 

Unfortunately, these experimental data are in direct conflict 
with published theory and experimental data for other 
viscoelastic materials. For example, experimental and theoretical 
solutions for simple isotropic viscoelastic solids allow for the 
existence of two creep compliance or relaxation modulus 
functions. As a consequence, the above approximations are 
application and material specific, and of limited interest. 

In the last 100 odd years closed-form solutions to a host of 
one, two, and three dimensional engineering design and scientific 
problems have been obtained using the theory of linear 
elasticity. Two and three dimensional solutions incorporating 
viscoelastic constitutive behavior have been far more difficult 
to handle because of the necessity of dealing with the 
f * iJkl (t kl/ a kl ) or f ij k i(t kl/ a kl ) . The differences between the 
generalized mathematical formulations and material specific 
solutions inferred from experimental data have also led to 
complications and confusion. 

The purpose of this document is to clarify the physics of 
material deformation and present a method for developing the 
simplest possible solutions for three dimensional constitutive 
equations of viscoelastic materials. Methods of application will 




also be illustrated. 

The three dimensional solutions may always be expressed in 
terms of tensor valued coefficients and a single scalar valued 
material dependent function, where the constitutive function has 
5 two independent responses, one time dependent and one time 
independent. The solutions for these constitutive functions 
represent a continuum in material behavior that also includes the 
theory of linear elasticity. Under the appropriate limiting 
conditions the generalized viscoelastic constitutive equations 
10 will degenerate into the linear elastic functions. 

^ Obtaining linear viscoelastic equivalents for two and three 

*? dimensional solutions defined in terms of linear elastic solids 
?i is now a relatively straight- forward process. For example, all 
\U solutions developed for the linear elastic solids may be adapted 

if* 

15 W to obtain solutions for the linear viscoelastic materials, 

a. 

13 regardless of the degree of material symmetry, and at least for 

S. 3 

|U static loading conditions. Additionally, quantification of 

Q constitutive functions is based upon experimental data obtained 

for each material of interest. Understanding the correct form of 
20 the equations greatly simplifies this task, reducing the time 

involved and associated costs. 
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SUMMARY 

This document presents a method for solving arbitrary multi- 
dimensional scientific and engineering design problems requiring 
solutions for stress, strain and deformation, and which therefore 
5 demand the incorporation of a material constitutive equation into 
the mathematical solution, and wherein the constitutive equation 
for solids, which quantitatively defines the relationship between 
stress and strain, incorporates two coefficient matrices of 
independent tensor valued coefficients and a scalar valued 
10 constitutive function, and where the values of the tensor valued 
«2 coefficients and the mathematical form of the constitutive 

•0 

function is specific to any particular material under 
[« consideration and the magnitude of deformation. 

w 

VS BRIEF DESCRIPTION OF THE DRAWINGS: 

15 13 So that the manner in which the above recited features, 

|U advantages and objects of the present invention are attained and 

m 

p can be understood in detail, a more particular description of the 
invention, briefly summarized above, may be had by reference to 
the schematic illustrations which are illustrated in the appended 
20 drawings. 

It is to be noted, however, that the appended drawings are 
conceptual illustrations of this invention and are therefore not 
to be considered limiting of its scope, for the invention may 
admit to other equally effective embodiments. 
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IN THE DRAWINGS: 

Fig. l.l - Molecular model of the bonded lattice structure 
of a solid material in undeformed and deformed configurations, 
and the coupling between the various components of an applied 
5 stress, a, and a resulting strain, e. 

Fig. 1.2 - The impulse elastic stress-strain and stress- 
ultimate strain curves a) , and the strain energy associated with 
ultimate strain b) . 

Fig. 1.3 - Idealized hypothetical stress-ultimate strain 
10 response of a solid material subjected to both elastic and 

;& inelastic deformation. So long as the material remains intact the 
*2 constitutive equations are always characterized in terms of 

H scalar valued constitutive functions, regardless of the state of 

Ifi 

Ul strain. 

IS 

15 10 Fig. 1.4 - Coincident reference frames x 1 - and z,- a), and the 

p reversal in the direction of action of shear stresses acting 

|U about the x^Xj and z 2 " z 3 ax ^ al planes caused by a 180 degree 

ffi . 

13 rotation of the z 3 reference axis with respect to the x 3 

M 

reference axis b) . 
20 Fig. 1.5 - The strain history of a viscoelastic material is 

directly dependent upon the stress history a^t), and controlled 

by the material functions f*(t) or f(t). 

Fig. 1.6 - Any arbitrary stress function a(t) can be 

expressed in terms of some initial stress a(0+) applied at time 
25 t = 0+, and a host of incremental increases in stress, the So- 9 

applied at the various interval times, the St 9 . 




Fig. 2.1 - Piece-wise continuous linear approximations to 
non-linear behavior are often reasonable approximations provided 
the range in stress and/or strain is appropriately restricted a) . 
As with linear materials, the stress-strain curve defines the 
5 components of elastic strain energy b) . 

Fig. 2.2 - A differential element of material loaded with 
the nine surface tractions, the o i} . The three surface tractions 
on each face may the be resolved into a single resultant 
traction, the t f . 
10 Fig. 2.3 - As the volume of a differential element 

approaches zero, Ax # Ay # Az -> 0, each of the surface tractions, 

P the a r} , and the resultants of those surface traction, the t f/ 

M 

represent stress at a point. The a„ and the t-, may then be 

y i 

p resolved into a single stress vector, a s , the state of stress at 
15 *3 a point. 

13 Fig. 2.4 - An orthogonal reference frame may be oriented 

|U with respect to the stress state vector, a s , such that one axis 

in 

O is collinear with a s , and a s is normal to the plane defined by 
the other two axes. 
20 Fig. 2.5 - The stress vector, a* v operates coincident with 

the z 1 axis, and normal to the plane defined by the z 2 and z 3 
axes. This stress vector is the solution for the state of stress 
in the z f reference frame a) . The state of stress solution, a* v 
for the Z| reference frame may represent the stress condition in 
25 the x f frame, defined by the stress vectors a 1 and <x 2 b) . 

Fig. 2.6 - As a consequence of an applied a* 1f there is an 
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associated set of strains, the 6* j1 , in the z i and x f reference 
frames. 

Fig, 2.7 - The stress state vector may also be defined by 
the a i vectors of the x { reference frame, and the direction 
5 cosines a, 6, and 4>. 

Fig. 2.8 - The relative contributions of non-linear strain e 
associated with any given applied stress a a) may be broken into 
a linear, e L , and non-linear component, € NL . The relative 
contributions of a non-linear strain solution for stress a b) may 
10 be broken into a linear, a L , and non-linear component, a NL , for a 
1^ given magnitude of strain, e. 

Fig. 3.1 - While the delayed elasticity functions f*°(a,t) 
and f D (€,t) are always non-linear functions of time, the ultimate 

^ strain solution may be either linear or non-linear a) . The non- 

1 

15 * u linear elastic impulse response, € E , and the delayed elasticity 
13 response, € D , may always be approximated by linear functions at 
|U time t->«, over an appropriately limited stress-strain domain b) . 

if? 

p However, these approximations will not necessarily be valid for 
any time t. 

20 Fig. 3.2 - The non-linear impulse response function and the 

ultimate delayed elasticity response functions may be described 
in terms of the sum of a linear component and a non-linear 
component . 

Fig. 3.3 - The delayed elasticity response of linear 
25 viscoelastic materials show strains, strain rates, and ultimate 
strains a) that are always proportional to the applied load c) . 
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At time t = oo the stress-strain curve is always linear. One 
possible subset solution for the delayed elasticity response of 
non-linear viscoelastic materials b) is a function that is 
proportional to some function of stress. € D (a,t) then takes the 
form of e D (a,t) = f*°(t)*g(a), where g(a) is a scalar valued 
function of the state of stress. At time t = « the stress-strain 
curve is then always non-linear. 

Fig. 4.1 - The cubic crystal lattice system of sodium 
chloride, NaCl, and lead sulfide, PbS. Bonding is orthogonal and 
equivalent in all three axial directions, 
p Fig. 4.2 - In additional to the cubic lattice structure a), 

2 orthorhombic lattices b) also have orthogonal bonding 

y 

^ arrangements. The constitutive responses of the bonds in both 

a - 1 ... 

!ff lattice structures, the f*-(t), may vary in each axial direction. 

IS 

Fig. 4.3 - The 'S' electron orbital is spherical shaped. 
13 The 'P' orbital shell is dumbbell shaped, and there are 3 
iy possible 'P' orbital directions, all of which are orthogonal. 

in 

p Fig. 4.4 - Chemical bonding commonly occurs through 

I* 

interaction of 'S' orbitals. However, if the 'S' orbitals are 
full they may act to repel each other and form anti-bonding 
orbitals a). Bonding produces a shared orbital shell b) . 

Fig. 4.5 - 'P' orbitals are more complex that 'S' orbitals, 
and there are correspondingly more bonding and anti-bonding 
configurations for these orbitals. 'P' and 'S' orbitals will also 
combine to form hybrid bonding orbitals. 

Fig. 4.6 - Shear stress, a 12 , within a plane of bonding a) 



may also be expressed as the single normal stress, o'^i if the 
reference frame is rotated 45 degrees. The racking of the lattice 
structure b) puts into play the forces associated with the 
stretching of the bonds, f, and f 2 ; the forces between the 
bonding and non-bonding/ anti-bonding orbitals resulting from 
angular distortion, f 3 ; and the reactive forces between the 
positively charged nuclei, f 4 . These forces all act 
simultaneously to produce a coupled response. 

Fig. 4.7 - If a force, a' 2 , applied to an orthogonal 
lattice racks the lattice in all axial directions, then all 

3 reactive forces operating between the bonds and nuclei come into 
\ play simultaneously. 

* Fig. 4.8 - A hypothetical bonding arrangement where the 

4 orthogonal bonds, between the atoms A-B and C-D, are connected 

^ via an intermediate set of bonds. The labeled arrows indicate the 
3 reactive forces produced by the electromagnetic fields and 

y induced by deformation of the lattice. 

fa 

r *. 

3 Fig. 4.9 - Many metal and metal alloys are characterized by 

face centered cubic packing of the atoms. Gold, silver, copper, 
lead, and many iron alloys have this bonding arrangement. This 
bonding arrangement produces symmetry in material properties in 
three orthogonal axial directions, the condition of orthotropy. 

Fig. 4.10 - Graphite, pure carbon, is characterized by 
stacked layers of hexagonally bonded carbon atoms a) . The layers 
of covalently bonded carbon atoms are connected by weak van der 
Waals bonds that act normal to the plane of each layer b) . Each 
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layer has a hexagonal lattice structure of interconnected benzene 
rings. 

Fig. 4.11 - The benzene ring is a hexagonal structure of six 
carbon atoms bonded by a resonating double-single covalent 
5 carbon-carbon bond. The plane formed by the ring has three axes 
of symmetry, where these axes are spaced 60 degrees apart. This 
particular degree of symmetry yields isotropic constitutive 
behavior in the plane of the ring. 

Fig. 4.12 - The body centered cubic crystal lattice system 
10 of metallic iron. While all the bonds between the iron atoms are 
q identical, none is orthogonal. 

Fig. 4.13 - Woven material with a 90 degree weave to the 
5^ individual fibers. 

5 - Fig. 5.1 - Newtonian fluids have a linear constitutive 

10 

15 ^0 relationship, while the constitutive functions for shear thinning 
!3 and shear thickening fluids are non-linear. 

IU Fig. 5.2 - The pressure forces, P, and viscous stresses, 

in 

13 T^, acting in the x-direction on a differential reference volume 
through which fluid is flowing at velocity v x . 
20 Fig. 5.3 - The velocity profile of steady state fluid flow 

in a confining channel a) , as well as the stress equilibrium 
condition for the static hydrostatic stress condition b) . 

Fig. 5.4 - A linear creep compliance and subsequent 
relaxation response of a uncross-linked polymer subject to a 
25 static stress a) . Only the deviatoric response shows 

unconstrained deformation under loading, b) illustrates the 
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deviatoric constitutive behavior of a Bingham material, defined 
in terms of applied shear stress and rate of deformation or 
strain rate. 

Fig. 6.1 - A cross-sectional surface schematic of a beam a) 
5 and the forces, the and bending moments, the tt rji acting 

within the beam b) . 

Fig. 6.2 - The moments and shear forces acting about a 
differential element of a beam, a) and b) , as well as the total 
differential resultant load, p»dx acting about the mass centroid 
10 of the beam, O, c) . 

p Fig. 6.3 - Plane sections, defined by cross-sections A, B, 

.p and C a) remain plane during deformation b) . The resulting 
j.q geometries of deformation are then used to define the curvature 
I« of the beam c) . 
15 ^ Fig. 6.4 - Two line positions a) within a beam retain the 

^ same relative separation during deformation b) . From the relative 

* N 

l|J changes in lengths of these line segments the resulting strain 

m 

P solution may be quantified as function of the distance y above or 
below the neutral axis. 
20 Fig. 6.5 - The symmetric distribution of the longitudinal 

loading function P x a) , and the differential distribution of that 
force b) which produces the resultant longitudinal force and 
the bending moment c) . 

Fig. 6.6 - The geometric relationship between the 
25 displacement w a) and the slope of the neutral axis that occurs 
as a resulting of bending b) . 
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Fig. 6.7 - A loaded cantilever beam a), the location of 
reaction force R c needed to support the beam b) , the location and 
direction of action of the maximum internal bending moment 
and shearing force c) , as well as the distributed internal 
5 bending moment and shearing force d) . 

Fig. 6.8 - Woven material with a 90 degree weave to the 
individual fibers. 

Fig. 8.1 - An arbitrary time dependent stress <x(t) and an 
associated delayed elasticity strain f* D (t) function. 
10 Fig. 8.2 - The accuracy of any approximation for a 

q historical delayed elasticity response is directly related to the 
^| time intervals between the incremental stresses approximations 
l n (see Figure 8.3). 

j| Fig. 8.3 - The accuracy of an approximation for a stress 

15 *9 history function, a) vs b) , has a direct bearing on the accuracy 

P of the strain history solution (Figure 8.2). The better the 

'"4 

iU stress history approximation the more accurate the solution for 

m • 

p strain history. 

Fig. 8.4 - For real materials a static stress, S,a, produces 
20 a time independent elastic impulse strain, e* t and a time 
dependent delayed elasticity strain, e^. 

Fig. 8.5 - For linear materials the incrementally applied 
static stresses, 6 n o, produce linearly additive time independent 
elastic impulse strains, the € n E . The time dependent delayed 
25 elasticity strains, e n D , are also linearly additive. All 

components of strain add linearly to give the total strain 
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response as a function of time. 



S3 
=g 

in 

IS 



!U 

in 

13 
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DETAILED DESCRIPTION OF PREFERRED EMBODIMENT 

SECTION 1 

LINEAR 3-D VISCOELASTIC CONSTITUTIVE FUNCTIONS 

One dimensional constitutive functions for linear 
5 viscoelastic solid materials may be written in the form of, 

o(t) - l • f (t) e , e(t) - 1 • f*(t) a 



where '1' is the unit vector and a and € are the magnitudes of an 
^ applied static stress or strain respectively. The relaxation 

modulus function, f(t), and the creep compliance function, f*(t), 
10 are scalar valued time dependent functions. a(t) and e(t) are the 

^ resulting time dependent stress and strain functions. The most 

IP 

05 generalized three dimensional Cartesian expressions for equations 



S3 1.1, may be written as, 



IVi 



m ait)^ - f ijkl (t kl ) e kl , €(t) ij . - f* ijk i(t kl ) o kl 



15 and imply complete independence for each of the scalar valued 
f Ukl (t kl ) and f* ijkl {t kl ) functions. These functions are analogous 
in form to the three dimensional linear elastic constitutive 
equations, i.e, €-- = D ljk i # <7 kl / etc. 

For very small strain conditions the stress and strain 

20 functions of equations 1.2 may be defined in terms of either 

material or spatial coordinates, with little loss in accuracy. 

15 



For convenience, and conforming to convention, we choose spatial 
coordinates. 

1.1 Notation, Mechanical Coupling and Generalized 

3-D Constitutive Functions 

Limiting interest to small infinitesimal strains, as in the 
theory of linear elasticity, e } and a- } may be defined in reduced 
indicial notation as, 
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P 

W 

m 



N 

III 
ill 
13 



and 



'ii ' 



'13 ' 



- o 



'31 



22 



'23 



°3 

a* 



- o, 



1.3 



'21 



6l " 




e 2 


~ e 22 ' 


e 3 - 


e 33 






- 2 




e 5 - 


2 e 32 , 




- 2 




1.4 


- 2 


e 31 ' 


e 8 " 


" 2 e 23 , 


e 9 


- 2 


G 2 i 





Expanding the second of equations 1.2 into a matrix format using 
reduced indicial notation yields, 



e,(t) 








e 2 (t) 






jf* 22 (t 2 ) 


€ 9 (t) 






/* 92 (t 2 ) 



f* a9 (t 9 ) 








f* 29 (t 9 ) 




°2 


1.5 


^* 99 (t 9 ) 




o 9 





where the 81 creep compliance functions, the f* ij -(tj), are assumed 
15 to be completely independent and are functions of only when in 
time, the t j# each a } is applied. 
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Coupling 

The assumed complete independence in the constitutive 
functions of equations 1.2 and 1.5 precludes reduction of the 
number of constitutive equations from nine to six since there can 
5 be no obvious symmetry in the assumed independent f* i j(t J ) or 

fjj(tj). Likewise, there can be no corresponding symmetry in the 
solutions for the shear strains or shear stresses. Assumed 
independence in the constitutive functions therefore produces 
solutions for €j(t) and o-(t) that are incompatible with 
10 established stress and strain theories, and which independently 
p establish symmetry in both the shear stresses and in the shear 
,p strains. The assumed independence in the constitutive functions 



!pj also violates the requirement for continuity in both the strain 

y and displacements fields imposed by strain theory. As a 

i.y 

15 v ~ consequence, the compatibility conditions, which demand coupling 

s 

P; between strain responses, are likewise violated. The 

!U compatibility equations for small strain theory are, 

In 



dx 2 dx 2 2 dxjdxi dx 2 dx 2 dx^ 

d*e 33 + ye 2 2 _ 0 JPe 2JL 



1.6 



dx 2 dx 3 2 dx 2 dx 3 



and 
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^33 + 






^£ 3l 


dx x dx 2 


dx* 


dx x dx z 


dx 2 dx z 




3**1, 






dx x dx 3 


dx 2 > 


8x 2 dx 3 


dx 2 dx 1 










dx 2 dx 3 


dx* 


dx 1 dx 2 
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Figure 1.1 illustrates bonding between atoms in some body, 
crystal or molecule of an arbitrary solid. Bonding and other 
interatomic forces distribute the induced internal stress produce 
5 by an applied load. These stresses in turn induce the distortions 
& between atoms in the molecule or the crystal lattice. However, 
? ~ deformation in one given direction, although time dependent, must 
!!j be coupled to deformation in any other direction. So, there is an 

if! 

W incompatibility between physical reality and the assumption of 

10 iP independence in the f*jj(tj) or f--(tj). 

Q It is also important to note that there is also an 

t . i 

Hj incompatibility between physical reality and the assumption of 
p the existence of 'perfectly linear' constitutive behavior. Aside 
from the non-linearity in stress and strain functions associated 

15 with changes in internal lattice geometries, material 

constitutive behavior is fundamentally non-linear due to the non- 
linear nature of the electromagnetic forces acting between atoms, 
regardless of the material type or bonding condition. Contrary to 
what empirical data may imply, any linearity in material 

20 constitutive behavior, elastic or viscoelastic, can only be an 
approximation where linearity implicitly imposes a number of 

18 



conditions that restrict 'linear' behavior to the small strain 
condition. 

1.2 Reduction of the Function Matrix 

Each scalar valued constitutive function, representative of 
a 'real' viscoelastic solid material, is always comprised of two 
physically measurable, independent and linearly additive 
component responses, 



f*(t> - f* B + f* D (t) 

1 • o 

g - 

The elastic impulse response of this arbitrary creep compliance 
10 !!! function, f* E , is time independent and a constant for linear 

5^ materials. The second term, f* D (t), is the time dependent delayed 

tfl elasticity response that is, in general, non-linear in time. 

P Expanding one of the nine strain terms of equation 1.5, say 



|U ^(t), yields a sum of component strains, 
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e x (t) - f* lt (t x ) a x + f* 12 (t 2 ) a 2 + ... + f* 19 (t 9 ) a 9 



15 r 1.9 

- e^it^) + e 12 {t 2 ) + ... + e 19 (t 9 ) 



and likewise for the other eight € f (t). Obviously, 

e tl (t x ) - fruity a x 

1.10 

and so on. If each static a. } is set to unity the f* ij -(t) are 
equivalent to each component strain function, e.g., l # f* 11 (t 1 ) = 

19 



e l1 (t l ). Variation in the magnitude of stress just scales the 
magnitude of the strain components, the e f j(t). So, equation 1.10 
may also be written as, 



1.11 



etc. 
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The Unreal Materials 

Consider the 'unreal' materials, those abstract and non- 
existent materials that lack an elastic impulse response. The 
f*ij(t) are now defined solely in terms of a time dependent 
delayed elastic response, f* ij (t) = f* D .j(t). Set the static 
a j = 1, and apply all nine at the same initial time. If the 
functions of equation 1.5 are evaluated at some time t 1# the 
f* D 1 -j(t) become constant valued, i.e., 









^12 


z 13 . 


■ z 19 








Z 22 


z 2i ■ 










^92 


Z 93 




and likewise 


at some 


other 


time 












z( 3 • 










& 


z' 23 . 


7* 






7 1 




Z93 


7 1 


where the z v 


and Z 


t 

U 


are the values c 



1.12 



1.13 



20 



respectively. One constant can be expressed in terms of another 
constant multiplied by a scaling constant, so the relationship 
between f* D n (t) and f*° 12 (t) at time t,, for example, becomes, 

2 12 - f* - k 12 Z xl - *„ f* D ii<ti> 1 ^ 14 

5 Likewise for f*° n (t) and f* D 12 (t) at time t 2 , 

Z[ 2 " f* D 12 i^ m ^12 ^11 " f+ D 11 it 2 ) ^ ^ 



n As time t 1 -> t 2 the values of the f^-jft) change, but the forms 
]q of the functions do not change. The function f* D 12 (t) at both 
^ times t 1 and t 2 is always defined in terms of f^^Ct). So, as time 



10 [{; t t -> t 2 , then f* 0 ,,^) -> f* u u (t 2 ) , which demands that Z n -> Z' n 



^0 Likewise f*° 12 <ti> "> f*°i 2 (t 2 ) , which demands that Z 12 -> Z' 12 . The 
O coefficients k 12 and k' 12 are time independent, so as t 1 -> t 2 
IU equations 1.14 and 1.15 produce, 

in 

13 



15 where k 12 = k' 12 for all time. Continuing for the remaining 80 
f*° fj (t) produces f* D ij (t) = k^f *° u (t) . This solution not only 
satisfies the requirement of coupling in time between each 
component of strain, it also satisfies the defining requirement 
for linearity at very large time, call it t = «, where each 

20 component of strain becomes a simple time independent and 
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invariant function, a constant, that is scaled by both the values 
of the delayed elasticity function evaluated at t = » and the 
applied stresses. 

The Real Materials 

Equation 1.8 is the definition of the creep compliance 
function for all real materials. Expressing one of the strain 
components, say e^t), of equation 1.5 in terms of each of the 
independent f* ij (t) yields, 

€ t (t) - fluUx) 0 1 + f* 12 (t 2 ) 0 2 + ... + f* 19 (t 9 ) 0 9 



^5 - e^itj + e 12 (t 2 ) + ... + e 19 (t 9 ) 



1.17 



10 •■!! and from equation 1.8, 

e 1 (t) - [f* s lx + f* D 11 (t 1 )]a 1 + ... + [f** 19 + f* D 19 (t 9 )]a 9 



m 
m 



1.18 



'"-4 

|U which may be rewritten as, 

in 
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e 1 (t) - f * * X1 o x + . . . + f * E 19 o 9 

+ f* D 11 (t 1 )o 1 + ... + f* D 19 (t 9 )a 9 



1.19 



or 



€l (t) - f*\<o j + f* D li (t J )a i ort 
15 1.20 



For all e^t) this leads to, 



22 



e^t) - tf* B ij) oj + [f* D ij {t j ))o j 

1.21 

If the f* ij (t j ) are evaluated at time tj = t = 0+, these 
functions must degenerate into the impulse solution. The 
resulting coefficient matrix may always be written as 
5 f* E ij s= m ij *f* E , due to coupling through bonding, where the are 
constants and m n = 1, by choice of definition. f* E is a scalar 
valued elastic impulse function, a constant. The impulse 
component represents an independent response, so eguation 1.21 

Y2 may be rewritten as, 

<0 

in 

|y 

*3 Each of the delayed elasticity functions may be written as 

E 

13 € D , s (ts) = e,,(t) - € E ,-5 = f* D ,.;(t;) when the a, = 1. From previous 

J;.j 'J J 'J 'J 'J J J 

|U discussions pertaining to the ' unreal' materials we know that 

m 

p when the f* D |j(tj) are evaluated at equivalent but different time 
I* 

15 intervals (equations 1.14 - 1.16) the matrix of delayed 

elasticity functions will degenerate into a matrix of constants 
and a scalar value function of the form f* D ij (t) = n-ff^lt) , 
where n„ = 1, by choice of definition. Equation 1.22 then 
becomes , 



20 3 3 3 3 1.23 
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The times of application, the t jf of each o } have no effect of 
the form of f* D (t), and the n^- need not be equal to the m^. Any 
equivalence between the m^- and the n u is a special case 
solution. The same procedures may be employed to derive the 
5 solution for stress producing, 

o.lt) - [a Sj f* + b u f D (tj) ]e, 



Equation 1.23 may also be written in terms of a single set 
of matrix coefficients and a single scalar valued constitutive 
i** function. Letting the t ; = t, for convenience, the expression for 
10 €^(t) from equation 1.23 becomes, 

l i 

111 e 1 (t) - a x + [m l2 f * B + n 12 f * D ( t) ] o 2 

I'd 1.25 
IS +...+ [m 19 f* E + n 19 f* D (t)]o 9 

13 where m,, = n t1 = 1. Likewise for the other eight 6j(t). Factoring 

'■■4 

IU out f* E + f* D (t) yields, 

m 

k » {t) f^+ft'U) 1.26 



15 where k 1t = 1. The right side of equation 1.22 now becomes, 

[ *« (t ' if - (t) "> - ' ^ f ;.'.:AB*' t) °' i.27 



Note that the k|j(t) are time dependent unless the = n^. 

Note also that kjj(t) is really k^-ttj) where there are six 
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different times of application of each static stress. There are 
therefore six different values for f* D (tj) and the ^(tp matrix 
is not symmetric unless the six t. = t, i.e., the six stresses 
are applied at the same initial time t Q . For example n^f* 0 ^) 9* 
5 n 21 »f*°(t 2 ) unless t, = t 2 even though n 12 = n 21 because the two 
stresses a, and a 2 do not have the same times of initial 
application (the symmetry proof for the matrix coefficients is 
given in the next section). Clearly, f* D (t 1 ) * f* D (t 2 ) and so 
kwCt,) * k 21 (t 2 ). 

10 Let's now evaluate equations 1.2 6 and 1.27 at time tj = 0+, 

p and tj -> oo. At time t } = 0+, f* D (0+) -> 0, these equations yield 
[§ the impulse solution where the k n (t) -> m ir When time is large, 
!~ say t- -> », then f* D (t->oo) approaches an asymptotic elastic 
modulus. The k^t-r^oo) are again constants, yielding, 

fP 

U «i(t— ) - t k ±j ] lf* B * f* D (t~-)}0 j - Iky] f*(t-0°)Oj 

15 1 • 2 8 

!U 

m ' 

O which is identical in form to the linear elastic constitutive 
I* 

equations, but only at this time limit. At these two limiting 
times (or at any time) the values of the k^-ft) are not equal 
unless the = n^-. 

20 The Number of Independent Constitutive Equations 

From stress and strain theories we know that shear strains 
are symmetric and the shear stresses are symmetric if the 
stresses are not coupled. So, there are only six independent 
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strains and six independent stresses. In the theory of elasticity 
the number of independent constitutive equations reduces from 
nine to six because of this symmetry. Evaluating at time t = 0+, 
and taking into account symmetry in the ij and kl indices of the 
5 m ij|a and a ijkl coefficients, the number of independent coefficients 
reduces from 81 to 36. Likewise for the b ijkl and n ijkl . The delayed 
elasticity component of equation 1.23 then reduces to, 

e D x it) 
e D 2 (t) 

e D 6 (t) 

Q 

,p and so on. 

in 

5:3 

10 I 1.3 System Energy and Symmetry in the Matrix Coefficients 

In the theory of elasticity symmetry in the coefficient 
H matrix may be demonstrated if a strain energy function exists. 
W The basic assumptions involved are that strain is totally 

in 

Q recoverable and deformation is adiabatic and/or isothermal with 
15 complete reversibility in heat transfer. Let's temporarily 

require deformation to be completely recoverable and slow enough 
such that it is isothermal, where heat transfer is completely 
reversible. Under these constraints a strain energy function will 
exist for the viscoelastic solids. 
20 The second law of thermodynamics may be written as, 



'13 



f* D (t 2 )- o 2 



1.29 



^61 ^62 ^3 



66 



f* D UJ- o fi 
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6E ti - U + j [JO + Q]dt 1.30 

U represents the mechanical energy converted to recoverable 
elastic strain energy, while JQ represents mechanical energy 
converted to heat during deformation. Q represents heat loss and 
5 E ti is total internal system energy. The integral term represents 
the net change in the total internal thermal energy, which has 
been constrained to be zero. For a complete deformation cycle 
this equation now reduces to, 

P 

| " i pdt i.3i 

in 

10 % where P is the stress power term, the time derivative of the 

iy 

elastic mechanical energy function, i.e, the strain energy 

E 

function U. 

m 
in 

U Symmetry in the Matrix Coefficients 

Both components of the function f(t) determine the magnitude 
15 of strain at any time t. Let's evaluate f D (t) at a large finite 

time interval, call it t->«, where the ' ultimate' strain has been 
reached. At this time the ultimate strain is the sum of the 
linear impulse strain, and the delayed elastic strain f D (t->co) 
(Figure 1.2a). Equation 1.24 now becomes, 
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0,(1:—) - ia ij f B + b ±j f D (t-~)] 6j - [ a' i:l + e., 

1.32 

where a' rj = a^f* and b' rj = b jj «f D (t->«») . Likewise for equation 
1.23, 

€.,(£:-«) - [ my f* B + ji^ f* B (t-») ] o_,- [m' i; , + 

1.33 

- + tVt—) 

where m',^ = m^f* 6 and n',j = n ij «f* D (t->») . 

Equations 1.32 and 1.33 are linear, so strain energy or 

complimentary strain energy functions derived from these 
'r solutions will also be linear. The strain energy function U is 
=** the sum of impulse component of strain energy, U E , and the 

in 



10 W delayed elastic component of strain energy, U D , i.e., 



m 

m 



0 V ' J ' 1 1-34 

iw 

in 

P (Note that deformation also produces a kinetic energy term, call 
it KE, that is equal to zero for static equilibrium conditions). 
If only the ultimate strain is considered (Figure 1.2b), the 
15 maximum possible value for strain energy is the area under the 
stress-ultimate strain curve. Given that stress and strain are 
monotonically related functions and that equation 1.34 is linear 
and continuous we may demand the equivalence, 
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1* JJ 



de J dej dej de , 

and likewise for 1^ and U°. 

Let's evaluate equations 1.3 5 for the derivatives of the 
strains e 1 and € 2 . From equations 1.3 2 and 1.34 the first 
5 derivative with respect to € 1 is, 



a 2 ^ 



£r " + ^''w* i. 3 6 



^ The derivative with respect to e 2 is then, 

2 - a 12 f* + 2> l2 f *(t~) , „ 

I s * 3e 2 12 12 1.37 

111 

.IS 

!9 Inverting the order of differentiation yields, 

13 

io iu a^a* 2,21 f (t } 1.38 

in 

13 

5 Equating the results of equations 1.37 and 1.38 then yields, 



a tl f* + b 21 f(t-~) - a l2 f + b u f>{t~) i 



or 



(a 



21 



a l2 ) f* + (i? 21 - b 12 ) f D (t-») - 0 



1.40 
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f (t) is non-zero in value for time t > 0, so the coefficient 
differences must equal zero, and a 12 = a 21 and b 12 = b 21 . These 
equivalences are independent of f E and f D (t) and therefore hold 
true for any time t > t=0+. Continuing for all possible 
5 variations will show that the a^- = a^- and b i} = b- when i 5* j . A 
similar result may be obtained from the complimentary strain 
energy function yielding m^- = m } - and n^- = n jM when i 5* j . Each 
matrix is symmetric, with only 20 unknown coefficients (recall 
that = n n = 1 by choice of definition) • Because symmetry is 
10 independent of strain and the constitutive function, it must also 
be independent of stress, and hence independent of the magnitude 
'? of the recoverable strain energy as well as any change in strain 

!_ energy. 

jfl 

W As a material deforms the change in total internal energy is 

ta 

15 18 a function of both the change in the elastic strain energy, and 
13 the change in internal thermal energy. However, if the elastic 
|U energy term in equation 1.31 is evaluated for symmetry, but at 

in 

13 some different temperature, we shall again obtain the symmetry 
" solution for the coefficient matrices. The values of constitutive 
20 functions and the matrix coefficients may be functions of 

temperature (see section 1.5), but matrix symmetry is not a 
function of temperature nor of a change in temperature. Matrix 
symmetry is independent of the state of internal energy, so 
coefficient symmetry is independent of both the state of total 
25 internal energy, and changes in total internal energy. 

Matrix symmetry is also independent of the state of strain, 
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including non-recoverable strain. Non-recoverable strain is the 
consequence of a different physical process other than simple 
elastic deformation, e.g., the breaking and reforming of bonds, 
as in ductile plastic deformation of metals. However, so long as 
5 the material remains intact, the viscoelastic constitutive 

equations at any strain state condition will not change form. No 
matter the magnitude of unrecoverable strain (Figure 1.3), the 
linear viscoelastic constitutive equations will always be given 
in terms of scalar valued constitutive equations and symmetric 
10 matrices of coefficients. The values of the matrix coefficients 
and the form of the constitutive functions may change, but the 
^ basic form of the equations does not. 

In 

\'y 1.4 Material Symmetry Considerations 

13 The number of non-zero matrix coefficients associated with 

15 13 any degree of material anisotropy may be determined using the 

U methods employed in the theory of linear elasticity, and require 



Ui 

in 

3 the vector transformation equations, 



dZ{ dz, dz: 

e *« " a^J c " ' ■ ST k ^ °« i-4i 



where the partial derivatives are the direction cosines. The 
20 constitutive functions f* E , f E , f* D (t) and f D (t) are all scalar 

valued functions and are therefore unaffected by reference frame 
transformations . 



31 



One Plane of Symmetry (Aelotropic Materials) 

A plane within a body of material is considered symmetric 
if, for every pair of coordinate systems located at a point on 
that plane, which are mirror images of each other, the matrix 
5 coefficients have the same values at that same point. If two 

Cartesian reference frames are defined to be coincident with each 
other, this symmetry condition is equivalent to flipping one axis 
of one spatial reference frame 180 degrees with respect to 
another axis of the other spatial reference frame (Figure 1.4). 
10 Define two reference frames, and z i , where the Xj-Xj and 

^ Z!-z 2 axes are coincident, but where z 3 = - x 3 . Let stress and 
? |= strain of the x f system be a and e, while stress and strain 

s - 

it associated with the z, system are a* and e*. The defining axial 
^ transformations are x, = z v x 2 = z 2 , and x 3 = - z 3 . So, equation 
15 (8 1.41 yields the following direction cosines, 

□ 



!U x ' Xt 1 1-42 



m 

13 



and 



dz 1 dz x dz 2 dz 2 3z 3 dz 3 
dx 2 dx 2 dx x 5x 3 dx 1 dx 2 



1.43 



Substituting equations 1.42 and 1.4 3 into equation 1.41 yields 
20 c* i - } = a i - j , and e* i} = e i} , except as follows, 
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°32' °*31 - " °X3 

1 . 44 



e+ 23 - - e 32 , e* 31 



Letting t- } = t in equation 1.24 produces the following two 
equations for a t (t) , 



and 



0^(0 - [a^f* + b 1:j f D (t) ]e*j 
_ 1*46 

□ 

■9 

!!! Substitution of the appropriate a* i and c*,- transformations into 



^ equation 1.46 produces, 

ID 

m 

s o x (t) - [ a^ + a 12 e 2 + a 13 e 3 - a 14 e 4 - a 15 e 5 + a 16 e 6 ]f 5 

^ + [ + Jb 12 e 2 + Jb 13 e 3 - i> 14 e 4 - i> is e 5 + i> 16 e 6 ]f D (t) 

ill 

10 J5 as well as the untransf ormed expression of equation 1.45, 

■ s 



o x (t) - [ a xl e x + a 12 € 2 + a 13 € 3 + a 14 e 4 + a 15 G 5 + a 16 e 6 ] f J 

+ [ + jb 12 62 + ±> 13 e 3 + 2> 14 e 4 + jb 15 e 5 + jb 16 e 6 ]f D (t) 

Subtracting one equation from the other yields, 

0 - [-2 a 14 e 4 - 2 a 15 e< ] f * + [-2 b 14 e 5 - 2 £> 15 e 5 ] f D (t) 



1.48 



1.49 



for all arbitrary values of € 4 and € 5 . Evaluating at time t > 0, 
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when the strains are not zero, equation 1.49 is satisfied only if 
the a 14 = a 15 = 0, and b 14 = b 15 = 0. Similar results will be 
obtained for transformations of a 2 (t) , a 3 (t), and a 6 (t) . 

Continuing for the shear stresses yields the following 
5 equivalences for a 4 (t), 



- o 4 (t) - [ a 41 e x + a 42 €2 + a 43 e 3 - a 44 e 4 - a 45 e 5 + a 46 e 6 ] f E 

+ [ jb 41 e x + i> 42 e 2 + 2> 43 e 3 - 2> 44 e 4 - ±> 45 e 5 + b 45 e 6 ]f D (t) 

and the untransf ormed expression, 

o 4 (t) - [ a 41 € x + a 42 e 2 + a 43 e 3 + a 44 e 4 + a 45 e 5 + a 46 e 6 ]f E 
^ + [ + ±> 42 e 2 + Jb 43 e 3 + jb 44 e 4 + i> 45 e 5 + £ 46 e 6 ]f D (t) 



1.50 



1.51 



Adding these equations produces, 

W 

10 0 - t 2a 41 + 2a 42 e 2 + 2a 43 e 3 + 2a 46 e 6 ] f E 

10 IS 1.52 
? + [ 2i> 41 € x + 2jb 42 e 2 + 2jb 43 e 3 + 2jb 46 e 6 ] f D (t) 

13 
H 

|U Evaluating as before yields a 41 = a 42 = a 43 = a^ = 0, and 
13 b 41 = b 42 = b 43 = b^ = 0. Similar results will be obtained for the 
transformation of a 5 (t) . 

The above results now produce the following solution for the 
15 a fj - matrix, 
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a ll 


a 12 


a i3 


0 


0 


a !6 






3 23 


0 


0 


a 26 


S 31 




3 33 


0 


0 


a 36 


0 


0 


0 


a 44 


a 45 


0 


0 


0 


0 


a 54 


a 55 


0 


S 61 


^2 


a 63 


0 


0 


3 66 



and similarly for the m^. and n^- matrices. The number of 

independent coefficients in each matrix has been reduced. to 13, 
because the = a jif etc. Note from equation 1.27 that the 
5 number of non-zero k^-ft) coefficients is also 13. This result 

— and solutions for the number of non-zero matrix coefficients 

O 

^ matrix for the other degrees of material symmetry or isotropy are 
!f identical to those obtained for the linear elastic materials. 

in 

lO Isotropic Linear Viscoelastic Materials 
10 □ For isotropic materials the constitutive equations 

IU degenerate to a form very similar to the equations of isotropic 

in ... 

=3 linear elastic materials. Skipping the derivation, which is 
straight-forward, the constitutive equations for the time 
dependent strain field reduces to, 

eft)- K * {t) o 6 + G * (t) o'.. 
15 e V (t) 9 °** 6 ^ 2 °^ 1.54 

where the <r kk m 6 iJ are the hydrostatic stresses, and the c' i} are 
the deviatoric stresses. The bulk modulus K*(t) and shear modulus 
functions G*(t) are defined as, 
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K*U) - 3(1 + 2^ 2 ) t + I I 1 " 12 ! f* D (t) ] 

(1 + 2m x2 ) i # 55 



and 



d - ™i 2 > 1.56 



Only when the m^- » n fj are the bulk modulus function and the 
5 shear modulus function described by what appears to be the same 
scalar valued constitutive function, f*(t). Similar expressions 
p may be derived for the time dependent stress field of equation 

'2 1.24 as well. 

!== 

in 

1.5 Temperature Effects 

10 '9 Both the constitutive functions and the matrix coefficients 

P may be functions, of temperature. Letting the variable 0 represent 

!U temperature equation 1.23 may be written as, 

in 
13 

H €i (t, 8) - lm ij (B)]f* B (B) oj + [n ij (Q))f* D (t, 0) o j 



The variation in the m^- and n^- coefficients as functions of 
15 temperature is easy to understand if we recognize that the bond 
forces acting between atoms need not be geometrically symmetric 
or orthogonal, and the bonds need not be identical or between 
atoms of the same element. The same holds true for the other 
interatomic forces. As a consequence, temperature variations may 
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produce geometric changes in lattice structure as different bond 
types, or similar types of bonds (covalent, ionic, etc.) between 
atoms of different elements expand or contract with different 
differential increments. Additionally, the relative strengths of 
5 each bond type, or of bonds between atoms of different elements, 
need not have the same differential variations for some 
incremental temperature change. Again, the nature of the forces 
acting between atoms demands that the expansion or contraction of 
the bonds be fundamentally non-linear. The same holds true for 
10 the other interatomic forces. So, the m l} {0) and n rj (0) are, in 
$□ general, non-linear functions of temperature. Furthermore, as 
h k shall become apparent shortly, they can also be functions of time 

as the temperature varies. 
|~ Now, for a body of material a non-uniform temperature 

15 10 distribution would demand non-homogeneity in material properties. 
G Restricting the temperature in a volume of material to be 
!U uniform, the m n (0) and n n (0) coefficients are then functions of 

m 

p the same independent scalar valued variable 9. If variations in 
temperature are sufficiently small such that there are 
20 inconsequential geometric changes to bonding angles and 

negligible differences in the relative changes of the bonding and 
non-bonding forces operating between atoms, then it is clear that 
the two coefficient matrices reduce to, 
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and equation 1.57 now reduces to, 



e^t. 6) - [m ij )f* B (6) oj * [n i:j )f* D (t, 0) o j 



Temperature Induced Strains 

In an undeformed condition let the state of internal stress 
5 be zero. Deformation will produce a change in internal thermal 
energy. Similarly, any change in material temperature will 
produce a change in material volume, i.e., it contracts or 
q expands. Return to equation 1.26 and set the applied stresses / 
Iq <Xj, to zero, change the temperature of the material, and allow a 
10 p corresponding change in material volume. This leads to a 
constitutive function of the form, 

ID 
IB 

L 6^,68) - A^t^e) 

U 1 1.60 

"-4 

1U 

in 

13 where 60 = 6 - 0 O . 0 is some new temperature, and 9 0 is the 
I* 

reference temperature where the functions A f (t,0) = 0, by choice 
15 of definition. The time t defines the time interval from the time 

of change from the reference temperature, 6 Q . 

Thermal expansion or contraction produces a stretch or 

compression in the bonds between atoms, as when subjected to a 

mechanical stress. The functions A f (t,50) must therefore, in 
20 general, be characterized like the mechanical constitutive 

functions, i.e., 
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AjU.SO) - ^^(58) + A» i U,bQ) £ ^ 

1.61 



Once again restrict the change in temperature such that geometric 
changes in lattice structure and differences in relative changes 
of the bonding arid non-bonding forces are negligible. Each 
elastic thermal strain function, the A E 1 -(6^) / and each delayed 
thermal strain function, the Ti°-{t,80), are functions of the same 
variables, and must likewise be coupled. So equation 1.61 reduces 
to, 



10 «i(t f «6) - Aiit.M) - a i A E (6B) + ^ £ A D (t t iB) 

- 1.62 

H 

tn 

10 j| where the a i and fi i are constants, and a 1 = 0, = 1, by choice of 
W definition. Note that this function also includes a subset of 
*3 solutions where the time dependent component is a product of two 
IU non-linear functions, i.e., 

in 

- A D (t,6Q) - A D (t) • H(6Q) 

1 • OJ 



15 The number of non-zero a- and j3 i7 and any equivalences 

between the a { , and between the 0., are determined from material 
symmetry considerations. For the special case of isotropy the 
coefficients a 1 = a 2 = a 3 , 0, = j3 2 = j3 3 , and a 4 = a 5 = a 6 = 0 A = 0 5 = 

0 6 = °- 

20 The functions that define equations 1.62 and 1.63 are, in 
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general, non-linear functions of temperature. If they are 
approximately linear for some range in temperature change, even 
for very large time, then equation 1-63 reduces to the 
approximation , 



A E and A D (t) are now scalar valued coefficients of thermal 
expansion/contraction. For a large time interval, t->°°, the 
function A°(t->») is constant, and from equation 1.59 the 
g function 1.64 reduces to, 

'? 

u e^t^SO) - [a i A E + P i A D (t-co)]6B ~ A i ( t—) 66 

10 ^ 1.65 

in 

?y 
m 

m 

c 

13 Linear Thermo-Viscoelastic Constitutive Equations 

\U Equations 1.59 and 1.64 produce the generalized thermo- 



i 



|3 viscoelastic function, 



1* 66 



15 where the time variable t has been used for the temperature 

function. Apply a set of stresses, o- t to a body of material. If 
the resulting deformation does not cause any appreciable 
temperature change in the material, the effects of loading 
induced deformation and temperature may be uncoupled (an 
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approximation) . The effects of temperature change and loading are 
now separate and linearly additive, so the resulting thermo- 
viscoelastic constitutive equation becomes, 



* e^t, 66) - [m ij ]f^ B o j * [n ij ]f* D (t) o j+ [a £ A B + p t A D (x) } 66 

1 • 6/ 



provided the temperature change has a negligible effect on f* E 
and f* D (t). Evaluating equation 1.67 at sufficiently large time, 
t->« and t->«, then f* D (t->«>) and A(t->«) become constants. 
Equations 1.28 and 1.67 then yield, 



*3 

Iq e^t— , 66) - [k^ ( t-«) ] f * ( t-~) a j + ^.(1-00)60 



1.68 



10 ^ which is identical in form to the Duhamel-Neumann equations of 

IS 

W linear thermo-elasticity where the coefficients of thermal 

O expansion are taken to be a simple constants. In fact, dropping 

'-4 

!U the time dependent functions in equation 1.67 when they are 

in 

13 negligible will produce the Duhamel-Neumann equations for 
15 'perfect' elasticity, a non-existent material condition. Similar 
functions may be written for the stress equations. 



1.6 Linear Hereditary Integrals 

Strain in any 'perfect' linear elastic material is solely 
dependent upon the magnitude of the stress acting on the material 
20 at any point in time. Stress may have been applied in different 
increments and at different times, but because the constitutive 
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equation is time independent the magnitude of the strain will 
always be known at any specific time, regardless of when the 
stress increments were applied. Viscoelastic materials don't 
behave in this fashion. The resulting strain is a function of 
time, even if the applied stress is static, so the strain history 
is very much dependent upon the stress history. 

Consider the three dimensional linear constitutive equation 
for static loading given by, 



€i (t) - Im^f^oj + [n ij ]f^(t j )o j 



10 q The strain component € n (t) is defined by, 



In 

jjj e xl (t) - % f* s a x + n 11 f* D (t 1 )a 1 



1.70 



O Let an impulse stress a y = ^(0+) be applied at time t = 0+. At 

4 

I.U another time t = t' some incremental impulse stress Sa^t') is 

m 

p applied. Figure 1.5 illustrates the resulting stress and strain 
15 solutions as functions of time. The delayed elasticity function 
is only a function of time, so the strains associated with each 
applied stress are additive and the solution for total strain is, 

e n (t) - f* B [ 0^0 + ) + 6 o^t') ] 

1.71 

+ n xl [ a 1 (0*) f* D (t) + 6 a 1 (t / ) f* D (t - t') ] 

There are two parts to this function, the initial impulse 
20 stress af6(t), and a series of incremental impulse stress 
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functions, the 6o'f6(t - t'). The strain at some time t' is the 
sum of all the strains caused by all the loading steps (Figure 
1.6) that have occurred at the various time intervals t', i.e., 

t 

e lx (t) - m xx f* B 0^0 + ) + m xl f* E j do x (t i ) + 1 72 

0+ 

n u ^(0 + ) f* D (t) + n^f f* D {t - t f ) dt* 

0+ 

5 where the differential stress, da^t')/ may also be defined as, 

13 

*| Integrating the impulse integral term equation 1.72 becomes, 
III e^U) - m ll f»* o^t) + n xl 0,(0+) £* D (t) 

W t . 1.74 

f.0 + „„ f f**(t- t'ji^iill dt / 

10 »♦ dt 

S3 Separate equation 1.74 into the impulse response e E 11 (t), 

m 

■ ^ « J 'll<« - IB,, f^O^t) 

10 n 1.75 



and the delayed elastic response c D 11 (t), 

e^U) - n lt 0,(0+) /*'(t) + u n f f*°(t - t') d °* ( f 7) dt' 

J at 



dt' 1.76 



The delayed elasticity term may be recast through integration by 
parts as, 
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# 



- iJ 11 o 1 (0+)f* D (t) ♦ n ll [f* D {t - t>) o^t') ] C 

0+ 1.77 

0+ 



When the bracketed term is evaluated, the sum of the second term 
and the strain term (0+) *f * D (t) is zero. So, equation 1.77 

reduces to, 

s €^(t) - Ji ll a 1 (t)f«(0 + > - o^t') df * D{ d t t ~ t§) dt' l lQ 

^ However, f* D (t) has no impulse response, so f*°(0+) -> 0 and 

equation 1,78 becomes, 

s • 

IU 6ll(t) " aF dt 1.79 

is 

S3 Now define the variable s as s = t - t', the elapsed time 

10 lli between the initial load and the first incremental load. With a 

in 

p sequence of incremental differential loads the elapsed 

differential time between each incremental differential load is 
simply dt', so ds = - dt'. From this definition the differential 
expression in the integral may be rewritten as, 

d£+ D (t - t 9 ) m . d£* D (t - t') 
15 dt' " d(t - t') 1.80 



Equation 1.79 then becomes, 
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_ d i f.\ „ t „ ,wv df* D [t - t') 



1.81 



Combining equations 1.81 and 1.75 produces the final solution, 
e ll(t ) - ^ ♦ n n j\( t ') ^fjf 1>82 

There are five more terms in the solution for the strain 
5 component c^t). Letting t- = t for convenience, and repeating 
this exercise for the remaining components produces, 



similar expressions will be obtained for the other five e f (t) , 

\u 

IP yielding the final hereditary integral solution, 

10 S3 6l(t) " * fj^ "ye - 1-84 

In similar fashion an integral expression may be derived for 
the relaxation phase, i.e., 

o j( t> - a i} f°* )( t) ♦ ^ / .,<*') d ^< C .- t /) /} dt ' 1.85 



where f D (t) is the relaxation modulus function and tj = t. 
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Temperature Effects 

As with the delayed elasticity response function in equation 
1.84, a time dependent temperature function S$(r) and the 
corresponding linear delayed thermal strain function A d (t) will 
5 produce a thermal strain response that is out of phase with the 
linear time independent function, A E . Following the same 
procedures as outlined above we may derive a hereditary integral 
function for the time dependent temperature function 60(7), the 
elastic thermal stain function A E , the delayed thermal strain 
10 function A d (t) , the coefficients a, and (3-, and time t. Beginning 
s 3 with equation 1.67, but skipping the details of the derivation, 
"S the complete form of the resulting linear thermo-viscoelastic 
i~ hereditary integral equation will be given by, 

w 

in t 

m - ♦ n Sj J aj (t') df ; ( 7 t .>) /) dt ' 1.86 

13 , 

%j + a i A B b%(x) + Pi f 66 (tO dA ji (T " V di' 

iy d(T " 

sn 

15 13 If the time dependent functions f*°(t) and a d (t) are 

H 

inconsequential in defining material response, equation 1.86 
again reduces to the Duhamel -Neumann equation of thermo- 
elasticity. 



ei (t,a6) - m ij f* B a j (t) + a i A B bQ(z) 

1 • 8 / 
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Sequential Loading 

When each dynamic load is applied at some different initial 
time, the time variables in the hereditary integral functions 
need to be adjusted appropriately. In such a case it may readily 
5 be shown that the hereditary integral function of equation 1.81 
may be written in the form, 

where the dummy variable k sets the order of application in time, 
t., of each successively applied loading function, the a(T k ). The 

13 

10 C* variables t. and t' are defined as, 



hi * * 1.89 

IS 
fa 



I** 



j A similar expression may likewise be derived for the stresses, 



y 
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SECTION 2 

3-D NON-LINEAR ELASTIC CONSTITUTIVE FUNCTIONS 

For a linear elastic impulse response the one dimensional 
constitutive relationships are, 

a - f * e , e - f* B o 

5 2 • 1 

where f* E and f E are the elastic constitutive functions, simple 
proportionality constants. Linear constitutive behavior is just a 
n limiting case to the larger, non-linear set of material 
% behaviors. As mentioned previously, linear constitutive behavior 

10 ^ cannot be an 'exact' solution for any material, and 'linearity' 

If] 

W is strictly a small strain approximation. For a non-linear 

f"S elastic response the one dimensional constitutive equations are, 

t 

^ a - f{e) e - f**(o) 

in 

s 

I** 

where f E <€) and f* E (a) are generalized non-linear impulse 
15 functions, representing all possible solutions, linear and non- 
linear (Figure 2.1a), large and small strains. 

2.1 The Power Series Representation 

Non-linear constitutive functions may be described in terms 
of 'nth' ordered power series functions, 
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o - f B (e) - c x e + c 2 e 2 + . . . + c n e n - c x e + £ c n e n 

n-2 2.3 

and 

z 

€ - f**(o) - ^ a + 2> 2 a 2 + . . . ♦ b n o n - ^ a + £ 2> n o fl 

n-2 2.4 

where f E L U) and f* E L (a) are the linear components, and the higher 
5 order terms are the non-linear components, f E ML (e) and f* E NL (<J) • 

The c k are arbitrary constants. These two power series functions 
n lack a free constant term, a c or a b . However, if there is no 
^ applied strain there can be no associated stress, or if there is 

'!« no applied stress there can be no associated strain, and so on. 

HI 

10 % Therefore c Q = b G = 0, always. When the higher order non-linear 
^ terms are negligible these equations reduce to equations 2.1. 

I U Linear Approximation 

s 

U ? 

O Because equations 2.2 - 2.4 are non-linear, incremental 

changes in stress and strain are not proportional to incremental 
15 changes in strain or stress. However, these functions may be 

approximated by piece-wise continuous linear functions within an 
appropriately restricted domain. If stress, a, is constrained so 
that a < a., then from Figure 2.1 the strain solution becomes, 



e - f* B (o) ~ a o 
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where 'a' defines the approximated slope of the stress-strain 
curve between zero strain and and so on. 



Power Series Sign Convention 

The power series functions are scalar valued magnitude 
5 functions, not vector functions. Recalling the definition of a 
vector, equation 2.4 becomes, 



€ - 1 • f* B (o) 



where 1 is the unit vector and the power series function f* E (a) 

\J 

% s is a scalar valued magnitude term. Stress, a, is likewise a 
10 scalar term. This definition insures that the vector is a non- 

in 

W trivial second order tensor. 

IS 

IS Different sign conventions define compression and tension. 

13 Tension is usually designated positive, with compression 

fij negative. This convention is fairly standard, but arbitrary. Now, 

in 

15 |3 the power series expansions contain both even and odd powers of 
stress and strain. Applying a unit tensile stress, a = 1, 
equation 2.4 yields, 



e - f* B (o) - Jb x + b 2 + Jb 3 + . . . + b n 



Applying a unit compressive stress, a = -1, negative in sign, 
20 equation 2.4 now yields, 
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Obviously, the two strain solutions differ in magnitude. By 
inverting the sign convention the strain solutions given by 
equations 2.7 and 2.8 should also invert in sign but maintain the 
5 same magnitude. However,, changing sign conventions has caused 

material constitutive properties to change, a nonsensical result. 
To achieve a proper solution, the scalar valued constitutive 
functions f* E (a) and f E (e) must be evaluated using the magnitude 

1*3 of the appropriate stress or strain, where the sign convention is 

if! . . . 

10 ;~ handled by the definitions, 



f* B (- a) -> - f* B (o) , f s (-e) -> - f E (e) 



IS 

IS 



S3 Let's again limit deformation to the small strain condition 

|U of the linear constitutive theory, but now allow for the 

in 

|g fundamental non-linearity in constitutive behavior for even very 
15 ' small strains. We may even allow the behavior to be wildly non- 
linear for very small strains. However, how 'wildly' non-linear 
the constitutive behavior may be for small strains is a function 
of the forces operating between atoms, and for small very strains 
is more an abstract notion than a possibility. 
20 Now, the small strain condition also implies negligible 

changes in the internal lattice geometry of a body of material as 
it deforms. Obviously, very large strains would rack the internal 
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lattice structure, causing directional changes in the internal 
stress distributions. The relationships between the resulting 
strains would therefore also change. By limiting, interest to 
appropriately small strains we allow for the fundamental non- 
5 linearity in the bonding and non-bonding forces operating between 
atoms, but preclude any significant changes in the internal 
structural lattice geometries. Because changes in internal 
lattice geometries are constrained to be negligible, any internal 
stress redistributions are also negligible and the internal 
10 stress distributions retain their original geometric 

proportionalities. It then follows that because lattice 
: B geometries and internal stress distributions remain relatively 
H unchanged the relationships between the various resulting 

m 

UJ directional components of strain will remain relatively unchanged 

15 Ifl and likewise remain simply proportional. 

i 

Si 

IU2.2 3-D Non-linear Elastic Constitutive Functions 

n From stress theory we know that if there is no coupling in 

s " the stress field then the shear stresses will be symmetric. From 
strain theory we know the shear strains are always symmetric. 
20 Admitting the physically inadmissable possibility that each of 

the strain or stress components of a constitutive equation may be 
completely independent and described by a separate constitutive 
function leads to the following three dimensional form of 
equations 2.2, 
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13 



'*li<«*l> 



2-10 



There are now nine constitutive functions for each of the nine 
solutions for stress and strain, representing some form of a 
stress state or strain state function. Equations 2.10 may be 
written in an even more general format, where each of the 81 
possible components of stress or strain are assumed to be 
represented by a single independent constitutive function, i.e., 



10 



O 
*B 

a - 

I* 

m 
y 

m 

a 

m 
m 

n 



2.11 



where once again we may define stress and strain in terms of 
spatial coordinates because of the small strain constraints. 
Using reduced indicial notation and expanding the second of 
equations 2.11 yields, 



«1 






f**22<«2> 


f* B 13 (o 3 ) . 
f* B 2i (o 3 ) . 


. . f* B 29 (o 9 ) 


2.12 


e 9 










.. f* B 99 (o 9 ) 





The assumption of complete independence in the functions of 
15 equations 2.10 and 2.12 precludes any consideration of coupling 
between the strain response or of symmetry in the solutions for 
the shear strains, a condition completely incompatible with 
strain theory. The indices i and j of equation 2.12 therefore 
range from 1 to 9, not from 1 to 6. 
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Each of the. f^-faj) or f^-taj) functions now represent a 
potentially independent strain or stress component. From equation 
2.12 the strain component e n may be defined as, 



Expanding one of the nine strain vectors, say e 1# then yields, 
e x - f*^^) + f**i 2 ( a 2 ) + ... ♦ f* B 19 (o 9 ) 

2.14 

s«l and likewise for the other eight € if as well as the a y of 
p equation 2.11. 

in 

W 2.3 Reduction of the Function Matrix 



10 !'S The nine functions in each of the nine columns of equation 

1.3 2.12 are defined in terms of the same independent variable. For 
1U example, the nine € n , 

in 

13 



are defined in terms of a r The solutions for the strain 
15 functions, the f* 6 !,^), are as yet unknown. However, these 

functions must be coupled due to bonding, and they are functions 
of the same independent variable, cr t . 

When the strain functions of equation 2.15 are evaluated at, 
say a y9 they degenerate into constants, 
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where the z n are the strains associated with a,. When these 
functions are evaluated at some other value of stress, a' t , they 
degenerate into another set of constants, 



One constant may be expressed in terms of another constant 
multiplied by another constant, so z 21 for example, may be 
expressed in terms of the strain term z n as follows, 

10 

m 

10 13 Likewise for all seven other z lM , which may be expressed as 

"'4 

IU z n = k ii #f * n^ a i^ # The k ii are siin P le scaling constants. 
q At o'y the strain term z' 21 may be written as, 



and so on for the other seven z' n . Driving a t -> a' v the values 
15 of f*^^) and f* E u {a % ) change, but the functions do not change. 

The scaling coefficients k' 21 and k 21 are independent so when a 1 -> 
o' y then f* E 21 (a 1 ) -> f* E 21 (a / 1 ) and f^Ca,) -> f* E l1 (a / l ). Equations 
2.18 and 2.19 then produce, 
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2.20 



~* ^21 " ^21 



for all a,-. Following the same procedure for the remaining seven 
f * E ,-i< a i) produces similar results and equation 2.12 becomes, 



1 ^**12(°2> 



*21 



f** 29 (o 9 ) 



/** 99 <o 9 ) 



1 



2.21 



where f* E 1l (a 1 ) is now a scalar valued function. Again, the simple 
proportionality between responses is obviously an approximation. 
As mentioned earlier, this approximation is allowed by the lack 
of any appreciable change in internal structural geometries, a 
consequence of the small strain condition. Continuing for each of 
10 (9 the functions in the other eight columns in the matrix yields, 



in 

w 

in 



111 

ifi 
Q 



x 21 



x 22 



^23 



*93 



K 29 



2.22 



The t* E i-{Oj) functions are scalar valued, and control all strains 
in all directions for any particular stress function. 

The above solution implies that each constitutive function 
15 is independent, but associated with a particular type of stress, 
i.e., normal or shear, as well as a particular loading direction, 
where more that one function operates in any direction. Given the 
requirement for coupling through bonding, and given that loading 
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from any direction results in the deformation of the same set of 
bonds, equation 2.22 makes no physical sense. The f* E V} (o } ) 
functions are scalar valued, so orthogonal reference frame 
transformations, which are completely spatial in nature, have no 
5 effect on the form of the functions. Furthermore, the non-linear 
functions of equation 2.22 must contain, as a limiting case, all 
possible linear solutions, where the linear constitutive 
functions are characterized by a single scalar valued 
constitutive coefficient. So, equation 2.22 can not possibly be 
10 the final solution. Each function must represent solutions to the 
, 3 same function, where that solution varies depending upon the 
"? magnitude of the stress and the sequence of application of each 
]Z stress. 

*.r s 

s- 3 

5M 

19 • 

10 2.4 The Final Reduction 

15 Q Application of Normal Stresses 

iy Let's now evaluate the nine scalar valued constitutive 

•H 

□ functions of equation 2.22 where each stress is applied 

separately and alone. A single stress then defines the total 
state of stress, a deviatoric stress. First, set eight of the 
20 nine stresses equal to zero, except a % . Equation 2.22 then 
reduces to, 

l 

The shear strains must be symmetric so there can only be six 



2.23 
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independent values for the kj,, and equation 2.23 must reduce to, 



x 2l 



2.24 



Applying a z and then a 3 separately and alone produces similar 
strain solutions and equation 2.22 becomes, 



l 



1 

t 22 



1 

<2 3 



2.25 



111 



m 

io IS 

13 



where the sequence of application of the stresses may be 
arbitrary. The different k^- come about from material symmetry, 
or lack of symmetry, considerations. 

Vector Transformation Equations 

Material properties must be examined as a function of the 
variation in the orientation of a spatial reference frame. The 
vector transformation expressions for stress and strain are, 



dx k dx 2 



dz H dz* 



dx k dx x 



' Jti 



2.26 



where the two reference frames are x f and z x . For Cartesian 
15 reference frames the spatial derivatives are the direction 
cosines. 



Single Plane/Axis of Material Symmetry 

Following discussions in section 1, define two orthonormal 
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coordinate systems, x i and z s , where the x^Xj and Z 1~ Z 2 axes are 
coincident and define a plane, but where z 3 = - x 3 . Stress and 
strain in the K i system shall be a and e, while stress and strain 
in the z i system shall be a* and €*. 

Consider the stresses and strains in the '2-3' and '3-1' 
directions. Solving for the direction cosines yields, 

dx x ' dx 2 ' 3x 3 2.27 

dz x dz x dz 2 dz 2 _ 3z 3 _ dz 3 _ 

dx 2 " dx 2 " &q " dx 2 dx 1 dx 2 

Substituting these results into equation 2.26 yields a*^ = o Vl , 
*0 and €*-- = € i; , except as follows, 

I* 

in 



10 sy 2.^0 



None of the k u coefficients in equation 2.25 associated 
iy with the normal strains may be eliminated, unlike the k^- 
H associated with the shear strains. Applying only a. the 
expression for the shear strain £ 4 in the x s frame is, 

15 " *" f * B » (0l) 2.29 



In the z i frame the equivalent expression is, 

_ *« f * 2.30 
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Flipping the z 3 axis and substituting the appropriate stress- 
strain transformations into equation 2.30 yields, 



-e 4 - * 41 f** ll (o 1 ) 



2.31 



Adding equations 2.32 and 2.33 produces, 



o - 2 * 41 fr'uloj 



2.32 



10 



□ 

i 1 ^ 



IB 
in 



ii j 

in 
13 



The constitutive function is non-zero for any non-zero value of 
a 1# so k 14 = 0. The same transformation result is obtained for € 5 , 
yielding k 15 = k 14 = 0, Equation 2.24 then becomes, 



i 

*» 
o 
o 



2.33 



Applying o z and a 3 separately, the same symmetry condition 
produces similar results yielding, 



A 31 
0 

o 

*61 



^22 
^32 

0 

0 
^62 



^23 

0 
0 



f**ll<«l) 



2.34 
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Two/Three Planes of Symmetry 

In addition to the axial rotation of z 3 = - x 3 for a single 
axis of symmetry, the rotations for 2 axes of material symmetry 
now include z, = - x t , z 3 = - x 3 , with z 2 = x 2 . The direction 
cosines for these equivalences are, 



dz 1 
dz x 



- - i, 



dz 2 
dx 0 



dz x 
3x7 



dz 2 



- i, 



dz 2 



dz 3 



dz 3 

^2 



2.35 



P 



Equations 2.26 then yield the following transformations of c 6 for 
the separately applied stresses a,, a 2 and a 3 , 



«12* " ~ e i2' °*1 " °1» °*2 - °2» °*3 " °3 



2.36 



10 IS 



Following the same procedures yields k 61 = k^ = k^ = 0, and 
equation 2.34 reduces to, 



iv 31 

0 
0 
0 



32 

0 
0 
0 



^33 

0 
0 
0 



f* B 12 (o 2 ) 



f* B l 3 (<>3> 



2.37 



Normal Stresses and Complete Material Isotropy 

Now consider complete material isotropy, beginning with the 
15 results just obtained. Again, the two reference systems are 

defined to be coincident initially. Rotate x- the reference axes 
90 degrees about the z 3 axis, such that the x 1 axis coincides 
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with the z 2 axis, the z 2 axis with the -x, axis, and the x 3 axis 
with the z 3 axis. Rotating the x 1 -x 2 axes through some angle a 
about the x 3 axis produces the transformations from the Xj-x 2 
plane into the z t -z 2 P lane given by, 



z 3 - x 3 , z x - Xjcosa + x 2 sina, z 2 - -x x sina + x 2 cosa 

2*38 



From equation 2.36 the 90 degree rotation of axes yields, 

0*1 * ° 2 ' °*2 " °1' °*3 " °3 

2.39 

□ 

% y Applying a. in the x, axis direction the expression for strain in 

^ the x. direction is, 
m 1 

iy 

10 IS 2.40 

s. 

P 

ill In the z 2 direction the solution for strain as a function of a* 2 

in . 

€ 2 * - ir 22 f** 12 <0 2 O 

2.41 



Rotate the z 2 axis so that it is coincident with the x 1 axis, and 
15 demand equivalence in the resulting strains (the condition of 
isotropy) . Substituting the transformations given by equation 
2.39 into equation 2.41 yields, 
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«1 - k 22 f * S 12<°l> 



2.42 



Equating equations 2.40 and 2.42 produces, 



*2 2 f**»(«i> - ^xK) 



2.43 



and the solution for f * E 12 (a l ) , 



A 22 



2.44 



!3 

a- is a magnitude term, so the scalar function f* E 12 (a) scales in 

^0 

H terms of the scalar function f* E n (a). Multiplying the 

in 

y coefficients of the second column of the matrix of equation 2.37 

LB 

Ifl by the constant term l/k^ yields, 



-•3 

IU 

ifl 

10 13 

Isfc 



1 




1 




^21 


1 


^23 




*31 


^32/^22 


*33 




0 


0 


0 




0 


0 


0 




0 


0 


0 
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Demanding the same degree of symmetry with regard to the J^-Xj 
plane equation 2.45 reduces to, 
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1 

*» 

0 
0 
0 



l/*aa 
1 

^32/^22 
0 
0 
0 



l/*33 

^23/^33 

1 

0 
0 
0 



2.46 



Redesignating coefficients equation 2.46 becomes, 

f* B (o x ) 



_ where m-- 

=0 



^11 


J"l2 




^21 






/n 31 


m z2 


jn lx 


0 


0 


0 


0 


0 


0 


0 


0 


0 


1, 


and 


f* E 



= = m^ = 1, and f* t l1 (a) = f* fc (a) 



2.47 



P 
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Isotropy and Matrix Symmetry 

When o y is applied, the distortional strains are € 2 and € 3 . 
Consider distortion within the 1-2 plane. The z 2 axis is rotated 
so it coincides with the axis, and equivalence demanded in the 
resulting distortional strains. In the x 1 direction the 
expression for the distortional strain € 2 is, 



€ 2 - m 21 f* E (o 1 ) 



2.48 



In the z 2 direction the distortion strain as a function of a* 2 
is, 



€ i* " ^2 f* E (<> 2 *) 



2.49 
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Substituting the transformations of equation 2.26 into equation 
2.49 produces, 



Equations 2.50 and 2.48 must be equivalent and yield the solution 
m 21 = m 12 . Likewise m 13 = m^. Similar transformations may be obtain 
for the distort ional strains in the 1-3 and 2-3 planes, and = 
m 31 , and so on. So, the m^- = when j ^ i. For complete isotropy 
equation 2.47 now becomes, 



In 



10 



1 


™12 


W 12 




^2 


1 


m 12 




m i2 


^12 


1 




0 


0 


0 




0 


0 


0 




0 


0 


0 





f* B (o 1 ) 
f* B (o 3 ) 



2.51 



10 L All these normal (but deviatoric) strains are functions of a 

a *1 single non-linear, stress dependent, scalar valued constitutive 

iu 

\U function. 



Transformation Inversion 

Inverting these transformations, and demanding a lack of 
15 equivalence in the resulting axial and lateral strains, the form 
and value of the constitutive function does not change, because 
it is scalar valued. Any lack of equivalence between strain 
solutions demands a lack of equivalence in the diagonal 
coefficients, i.e., m n ^ ^ m^, and off-diagonal coefficients, 
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In fact, inverting all the previous axial symmetry 
transformations has no effect on the constitutive function (the 
physics of this result, as well as exceptions to it, are 
discussed in section 4) . Regardless of the degree of symmetry, 
the material is described by a single scalar valued constitutive 
function, at least with regard to the three normal strains, which 
are deviatoric because the a f are deviatoric, having been applied 
separately and alone. For complete material anisotropy equation 
2.25 becomes, 



10 □ 



15 



20 



%i ^12 

™21 ^2 m 22 

m 21 /n 32 m. 



33 



f**(o x ) 
f**(o 2 ) 



f**(o 3 ) 



2.52 
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m 



™€1 ^2 ^3 

Matrix symmetry for any degree of material symmetry is yet to be 
proven . 

Shear Stress and a Single Plane/Axis of Material Symmetry 

From stress theory we know that any shear (deviatoric) 
stress or shear strain vector may be converted into a normal 
stress or normal deviatoric strain vector through the proper 
change in the orientation of the spatial reference frame. 
However, before making use of this property, let's determine the 
number of non-zero matrix coefficients associated with a shear 
stress. Set eight of the nine stresses to zero, except a 4 . 
Equation 2.22 then becomes, 
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2.53 



13 

in 



where there are only six independent k j4 due to symmetry in the 
shear strains. 

Define two orthogonal reference frames, the x f and z f , to be 
initially coincident. For the single axis of symmetry the axial 
transformations are x, = z,, x 2 = z 2 , and -Xj = z 3 . Let <r f and e i be 
associated with the x { frame, while <x*j and e* i are associated 
with the Zj frame. The solutions for the direction cosines are 
given by equation 2.27. The stress and strain transformations are 
a., = and e rt = €*,-j, except for those of equation 2.28. 

From equation 2.53 the normal strain 6 t in the x f frame is, 



m 
m 



e 1 - k lt f** 14 (o 4 ) 



2.54 



y 

'4 



IU In the Zj frame the equivalent expression is, 



in 

13 



*1« '14(0*4) 
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15 Substituting into equation 2.54 the appropriate transformations 
produces , 



e 1 - -k lt f* B lt (o t ) 



2.56 



Subtracting equation 2.54 from equation 2.56 yields, 
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2 * M f*\<(o A ) 



2.57 



f*(a 4 ) is not equal to zero for any non-zero value of a 4 , so 
k u = 0. Similar results will be obtained for the transformations 
of the other two normal strains and k 24 = k^ = 0. No information 
is obtained for the coefficients k u and k^, other than they are 
not equal to zero. However, this transformation does yield 
k^ = 0, and equation 2.53 reduces to, 



10 



c - 

I* 

in 

w 
is 

19 
t 

13 
M 
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fa 

m 



0 
0 
0 

*54 



2.58 
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Following the same procedures for the other shear stresses and 
resulting strains will produce similar results. 

Shear Stress and Two /Three Planes of Material Symmetry 

The axial transformations for two/three axes of material 
symmetry are z 2 = x 2 , z 3 = -x 3 , and z 1 = -x v Solving for the 
direction cosines yields equation 2.35. This symmetry condition 
yields no new information pertaining to k u . However, for € 5 in 
the X| frame we have, 



e 5 - k S t f** 14 (o 4 ) 



2.59 



and in the z s frame, 
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Substituting the appropriate transformations for stress and 
strain into equation 2.60 yields, 



«s - -* 54 f.< 14 K) 2m61 



5 where, by definition, f* E (-a) = -f* E (a). Subtracting equation 
2.59 from equation 2.61 produces, 
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1:9 
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-2 ic S4 f* K 14 (a A ) 



2.62 



f*(a 4 ) is not equal to zero for any non-zero value of a 4 , so 
Jc^ = 0 and equation 2.58 reduces to, 



o 

0 
0 

0 
0 



2.63 
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The other shear stresses, a 5 -> a 9 , yield similar results for the 
other sets of coefficients. 

Now, rotating a reference frame 45 degrees will convert any 
shear stress into a normal deviatoric stress, although the 
magnitude will change. Given that f* E u (..) is scalar valued the 
function is unaltered by the transformation (obviously, a spatial 
reference frame transformation does nothing to alter physical 
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material properties) . A normal stress produces a normal strain, 
so f* E u (..) must be equivalent to f* E (.,)* The same equivalence 
must also hold true for the remaining f* E 15 (..) -> f* E 19 (.0* 
When the £ i are evaluated for a 4 , or any other shear stress, no 
5 other non-zero values for the X rj are found for the condition of 
complete material isotropy. However, we do find that 

Complete Material Isotropy 

Symmetry in the shear stresses yields equivalence in the 
10 solution of the constitutive function for each of these stresses, 

t —a 

^0 i.e., f* E (a ii ) = f* E (a H ). When the are evaluated six shear 

H stresses yield no more information about the strain field than 

!U three shear stresses, so equation 2.12 reduces to, 



€ 1 




1 


^12 


m X2 


0 


0 


0 






«2 




^12 


1 




0 


0 


0 




f**lo 2 ) 


«j" 






^12 


1 


0 


0 


0 






€4 




0 


0 


0 


^44 


0 


0 




f**(a 4 ) 


€5 




0 


0 


0 


0 




0 




f**<a 5 ) 


*6 




0 


0 


0 


0 


0 


/n 44 




f**(a 6 ) 



2.64 



15 where the stresses are applied separately and alone (recall that 
for every operating there is also a a- operating and that 
a v s a }i ) . The coefficient matrix is identical to that of linear 
elastic materials. 

As discussed, we may invert all the previous reference frame 

20 transformations and material symmetry constraints, and arrive at 
a 9 x 9 matrix of constant coefficients. However, there are only 
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six equations/ six unknowns, and 36 unknown coefficients, so the 
most general form of the equation is, 
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As a general condition matrix symmetry has not yet been 
5 demonstrated. However, the matrix coefficients must be symmetric 
in order to satisfy the limiting condition of linear elasticity. 

■y Linear Approximations 

H Non-linear constitutive functions may be approximated by 
\U linear functions, where the logical choice to begin the 
10 is approximation is the unstressed-unstrained condition. The non- 
13 linear constitutive function will be determined from the 

''4 , . 

IU application of one applied stress function, say a normal stress, 

\h 

i«* which completely defines the state of stress. Using such an 

approximation the solutions for the strains induced by different 
15 loads may then be superimposed, as with the linear elastic 

materials. However, once the total state of stress, the resultant 
of all stresses, exceeds the point on the stress-strain curve 
where a linear approximation is no longer valid, the resulting 
strain solutions become increasingly inaccurate. 
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2.5 Multiple Loads and the Non-Linear Solution 

The Stress Field and the State of Stress 

Recall that there are nine possible stresses, the a^-, that 
may act about some element of solid material (Figure 2.2). Those 
5 nine vectors may be expressed in terms of just three resultant 

stress vectors, t t , t 2 , and t 3 , each of which is the resultant of 
the three stresses acting on one face of the element. Shrinking 
the volume of the element to where Ax*Ay»Az -> 0, then the area 
of each surface face also approaches zero and yields the solution 
10 for stress at a point, a s (Figure 2.3). This resultant stress has 
a magnitude and a direction of action with respect to the primary 

13 

*0 reference frame which may be calculated. 

*D 

H Again orient an orthogonal reference frame with respect to 

ill 

UJ the vector a s , such that one axis, say the z, axis, is coincident 

ID 

15 IQ with a s , and the z 2 -z 3 plane is normal to a s (Figure 2.4). The z i 
§g system is now a reference frame for which the maximum state of 
sij stress is known, i.e., the type of stress, the magnitude, and the 

in 

£«3 direction of action. 

M 

Also again recall that there is also a set of orthogonal 
20 axes, oriented in some fashion to a primary reference frame, 

where the nine stresses may be expressed in terms of just three 
orthogonal normal stresses. These three normal stresses, the 
'principle normal stresses', are coincident with the axes of this 
new reference frame, the 'principle axes', and are just another 
25 representation for the state of stress. 

72 



\ 



Matrix Coefficients and Reference Frame Transformations 

The values of the m iikl are, in general, dependent upon the 
orientation of the spatial reference frame with respect to a body 
of material. The m ijkl matrix represents a fourth order tensor and 
the values of each m ij - kl coefficient vary according the 
orientation of respective reference frames, transforming 
according to the function, 



dx a dx b dx k dx 2 

ljkl " a^a^ a^a^ m&bcd 2.66 



where the partial derivatives are the direction cosines in 

□ 

S Cartesian coordinates. If the m iikl are known for one choice of 

0 

^ reference frame, then the values of the m ijkl in another reference 

n 

ij frame may be determined from equation 2.66. For complete isotropy 

i 

g the m* ijkl = m ijkl for all reference frame orientations. 

i 

HI Isotropy and Non-Linear Elasticity 

Load a differential element within a solid body such that 



the only stress a* 1 is oriented along the z t axis of the 
orthogonal z i reference frame (Figure 2.5a). a* t represents the 
solution for total stress, both in magnitude and direction, and 
therefore represents a s . Orient an x f reference frame initially 
coincident with the z % reference frame, then rotate the x f system 
about the x 3 axis (which remains coincident with z 3 ) so that 
z 1 -z 2 axes are oriented with respect to the x^Xj axes by some 
angle of rotation a (Figure 2.5b). In the x^Xj plane the stress 




vector a* 1 may be equivalently expressed as the resultant of two 
simultaneously applied component normal stresses, a t and cr 2 , 
operating coincident with the x t and x 2 axes, a*, now represents 
the solution for the total state stress, a s , in the z % system. 

Isotropy insures that the m ij|a are invariant for any 
reference frame orientation. Now, the generalized governing 
constitutive equation of non-linear elasticity is 
€ s = [nijjl^f * E (aj) . For the condition of complete isotropy the 
solution for the constitutive equation is, 
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where the individual component stresses are constrained to be 
applied separately and alone. 

From equation 2.67 the loading function a*, produces only 
three strains in the z f system (Figure 2.6a), €*, = e* 11# 
e* 2 = €* 21 and €* 3 = c* 31 . The solutions for € 1 and e 2 , associated 
with the and x 2 axes, are obtained by adding the appropriate 
vector components of €* n and €* 21 (Figure 2.6b), 



2.68 



^ 2 - e* 11 sma + e* 21 cosa 



From equation 2.67 the solutions for the strains, €* 1 and €* 2 in 
20 the z f system are, 
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e *l " " ^11 f * B 

2.69 

Substituting equations 2.69 into equations 2.68 produces the 
solutions for c 1 and e 2 in terms of the state of stress, 
a s = a* v 

€ i " cosa f**(o* t ) + m 12 sina f**^^) 

5 2.70 

€ 2 - /n lx sina f**^*!) + /n 12 cosa f**^*^ 

So, only when f* E (a s ) may be approximated by a linear 
function may the base reference frame solutions for strain be 
*Q added in a linear fashion to obtain a complete solution. A linear 

sO 

H approximation of f* E (a s ) will yield reasonable solutions for 

in 

10 Id strain from the component stresses only if the solution for the 

19 

ijg state of stress does not exceed the range of accuracy of the 

£ 

13 linear approximation (Figure 2.1a). 

m 

in 

]nj Generalized Non-Linear Solutions 

To obtain explicit strain solutions for simultaneously 
15 applied loads first define the stress components, the a } , with 
respect to a base Cartesian reference frame, call it x f . Next, 
solve for the magnitude and direction of action of a s with 
respect to the x f system. Orient another Cartesian reference 
frame, the frame of state z., such that one axis is coincident 
20 with a s , and a s is normal to the plane defined by the other two 

orthogonal axes (Figure 2.7). Next, define the direction cosines 
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for the Z| system relative to the system, and compute the 
transformation solution for the matrix coefficients m^-, as given 
by equation 2.66. Finally, solve for the strains in the z f system 
with a constitutive function of the form, 



where the index j may be either 1, 2, or 3, and a s is a normal 
stress vector aligned coincident with one of the z^ axes. The 
index i will vary from 1 to 6. When solutions for the €* f have 
been obtained in the z- frame they may then be transformed into 

□ 

10 *8 solutions for strain in x_- frame with the vector transformation 

=0 

H equations 2.27. 

sn 

: : 

" 5 i 

i 

ifl An Alternate Form of the Non-linear Solution 

13 In the z i frame f* E (a s ) may always be written as a power 

iy series expansion in a s , 

j* f*Ho a ) - a x c g + a 2 a/ + . . . + a a a s ° - £ a R a g » ■ 
15 n-i 2.72 

and may always be expressed as the sum of linear component, 
f* E L »a 8 , and a non-linear component, f* E ML <<* s > • The solutions for 
the e* i may then be written as, 
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e*, - [jn* il ]/*»(a.) - [m* n ] [f* B L o s ♦ ^'^(o,)] 

2.73 

where the resulting strain field has two components, one linear 
and one non-linear. f* E (a s ) reduces to the first order term if 
the material is linearly elastic (an approximation) . f * E L then 
5 becomes f* E . 

2.6 Symmetry in the Coefficient Matrix 

It was possible to show that for linear materials the m ;] 
are symmetric in the indices, where the m^- = m^, when i ^ j. 

P 

*S This symmetry produces a maximum of 21 matrix coefficients for 

*B 

10 M anisotropic materials, where m., = 1, by choice of definition. 

in 

fij Coefficient symmetry for complete anisotropy also established 

18 

|p symmetry for all other possible degrees of symmetry and isotropy, 
p valid for all possible choices of reference frames. 

m 

The Strain Energy Functions 

15 ^ Figure 2.1b illustrates the one-dimensional form of the 

stress-strain relationship of equation 2.73. For any value of 
stress a, there are two components to the solution for strain €, 
e L and € ul (Figure 2.8a). Similarly, for any value of stress e, 
there are two components to the solution for stress a, a L and 

20 (Figure 2 .8b) . 

For either linear or non-linear elastic materials the work 
required to get from the undeformed state to a deformed state of 
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a unit volume may be expressed by the one dimensional integral 
expression, 



W - U - J o de 



2.74 



W is work and U is strain energy. (Note that the work function 
5 defined here assumes static equilibrium, where the kinetic energy 
component, KE, is equal to zero. During the dynamic process of 
deformation the work function is given by W = U + KE.) The 
complementary strain energy function, U*, is defined by, 

:D J 2.75 

in 

10 ii] In three dimensions equation 2.77 becomes, 



fg 

1.3 



U* - f « y d - / €j d o. 



15 



Q -o 2.76 

!L; Now rewrite equation 2.73 in a one dimensional format and 

^ substitute the result into equation 2.76. This produces the 
solution for complimentary strain energy, 

U* - I e do - | (e^ + e^do 2 77 



which is now the sum of linear and non-linear components. 

Let's now expand equation 2.77 into a three dimensional 
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format by substituting the generalized form of equation 2.73 into 
equation 2.76 , 
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This function expands into, 



U* 



In the z { frame there is only one stress function, o 8 , which is 

P 

*B equivalent to a v so equation 2.79 reduces to, 



!? 

ID 

10 

i 

= n 

Sea 

iy The Linear Energy Component 

m 

10 ;^ Recall that in the z f frame the matrix coefficients m* ij|a are 

^ transformations of the m ijkl from the base reference frame 

(equation 2.66). If that transformation is inverted we must again 
obtain the same m ijkl . Furthermore, the product m* jJkl *f* E L will 
transform into m ij|cl *f* E L because f* E L is scalar valued and 

15 unaffected by coordinate transformations. An equivalent solution 
for the linear component of the strain field in the base 
reference frame may therefore be written as, 
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where the a- } represent the stress components in the base 
reference frame. This equivalence is obvious when we recall that 
when the non-linear term in the constitutive function f* E (a s ) is 
5 negligible or non-existent (which is never the case) , there is 
only a linear solution for the state of strain. 

The linear component of the complementary strain energy 
function may now be defined in terms of the. m f - kl , f* E L , and the 
_ component stresses , Oy Because the strain field must be 

10 *S continuous and non-zero for any given applied stress field, the 

=0 

^ linear component of the complimentary strain energy function must 

in 

W exist for all possible non-linear materials (i.e., every 

ifi 

|fl material) . The following equivalence must therefore hold true for 

i 

13 the linear component of the complimentary strain energy equation, 



15 M de^j de J de 1 2.82 

H 

Skipping the algebra we find that the = m.., when j 9* i, for 
all degrees of material symmetry and isotropy. The solution for 
symmetry is also independent of strain, stress, the constitutive 
function, and system energy. 

20 2.7 The Non-Linear Constitutive Equations for Stress 

The constitutive functions for stress may likewise be 
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expressed in terms of a matrix of coefficients and a single 
scalar valued constitutive function, 



a 22 
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where the a fj are symmetric in the indices ij, and a„ = 1, by 
definition. The number of non-zero a fj for any degree of material 
symmetry and isotropy are the same as for the m^-. When written 
in the above format, the individual strains in equation 2.83 are 
constrained to be induced separately and alone. 

The method of solving for the resulting stress field is as 
outlined for the strain field, except that the constitutive 
function f E (e) is a function of the strain state vector, e s , 
oriented in some direction to the primary reference frame. The 
solution for equation 2.83 in the frame of state then becomes, 



- [a* ii ]f*(«.) 



2.84 



where the a*^- have been transformed from the base reference 
frame to the frame of state. Once the a* i have been determined in 
the frame of state, the a f in the base reference frame may be 
obtained from the vector transformation equation 2.27. 
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SECTION 3 

NON-LINEAR 3-D VISCQELASTIC CONSTITUTIVE FUNCTIONS 

From the previous section it follows that a one dimensional 
constitutive function for a non-linear viscoelastic material may 
5 be written as, 



e(a,t) - l-f*(o,t) o(€,t) - l-f(e,t) 

j « J. 



where 1 is the unit vector. f*(a,t) and f(e,t) are scalar valued 
non-linear creep compliance and relaxation modulus functions 

§*"! 

respectively, and are functions of time and a static stress, a, 

10 H or a static strain, €. Like linear viscoelastic materials the 

=n 

i t 5 

uj non-linear viscoelastic materials must be described as the linear 

SB 

10 sum of two functions that describe separate, measurable, and 

q physically distinct responses - a time independent elastic 

jij impulse response, and a time dependent delayed elastic response. 

in 

15 j«* This should be obvious because the forces operating between 

^ bonded atoms in a body of material are fundamentally non-linear 
and the linear viscoelastic functions can therefore only be 
approximations. These two component functions add linearly to 
give the complete solution at some time t > 0, i.e., 



f*(o,t) - f* B (o) + f* D (o,t) 



and 



82 



f(e. t) - f(e) + f°(e.t) 



Applying the assumptions, restrictions and methods of the 
previous section let's now discuss non-linear viscoelasticity . 

3.1 The Non-linear Constitutive Equations 

5 Allowing complete independence between each possible 

constitutive response, equations 3.1 and 3.2 may be expanded into 
the three dimensional constitutive equations, 

e^a^, t) - f* ijkl (o kl . t kl ) - f* B ijkl (o kl ) + f* D ijk2 (o k i* t kJ ) 
*0 3.4 

in 

\il and 
(0 

10 p 3.5 

IU 

€ E ijkl (^ k i) and c° ijkl (a kl# t kl ) are the impulse and delayed elastic 
^ strain components respectively. Likewise for the stresses. The 
shear strains are symmetric, as are the shear stresses unless 
there is stress coupling. As mentioned previously however, the 
15 assumption of complete independence in the component functions is 
incompatible with the requirement for symmetry in both the shear 
strains and the shear stresses, and incompatible with the 
compatibility conditions of strain theory, which demand coupling. 
Using reduced indicial notation the expressions for each 
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in 



i!5 

5 ^r 1 



component of equations 3.4 and 3.5 e^la^t) for example, may be 
written as, 



e 1 (a jt t) - e^t.. ) + e 12 (..) + ... +«„(..) 

3.6 



where € n ( . . ) is, 
and so on. 



3.2 Sign Conventions 

As discussed previously, the non-linear impulse response may 
be defined in terms of a power series expansion. Power series 
10 W expansions contain both even and odd powers of stress or strain. 



M A standard but arbitrary sign convention is to define tension as 



positive, while compression is negative. To achieve a proper 



B solution to these functions the f* E (a) and f E (e) must be 
evaluated using the magnitude of the appropriate stress or 
15 strain, where the sign convention is handled by the definitions, 



f* E (- o) -> - f**(o) f B (-e) -> - f B (e) . 

3 • o 



If the delayed elasticity component of the non-linear 
viscoelastic constitutive function is evaluated at very large 
time, designated t -> «, the ultimate strain will have been 
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reached (Figure 3.1a). The one dimensional stress-ultimate strain 
function then takes the form of, 

e(0,t-«) - f**(o) + f* D (o,t—) 

Now, €(a,t->») may be expressed in terms of a power series 
expansion in stress, like the impulse strain component f* E (a) . 
So, f* D (a,t->«) may also be expressed in terms of a power series 
expansion in stress. Given that there may be terms of both odd 
and even powers of stress in f*°(a,t) that produce the solution 
f * D (a, t->«) , we likewise impose the sign convention, 

f* D (-o,t) -> -f* D (o,t) f D (-e,t) -> - f D (e,t) 

3.3 Partial Reduction of the Delayed Elasticity Function 

Matrix 

The Strain Solution 

From section 2 we known that for the small strain condition 
the non-linear impulse response may be expressed in terms of a 
scalar valued constitutive function and a matrix of coefficients. 
This component of strain is time independent, so equation 3.4 may 
be rewritten as, 
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where each stress is now constrained to be constant, and to be 
applied separately and alone. 

Expanding the delayed elasticity term in equation 3.11 in 
matrix format yields, 



e D i:f (o jt tj) - 



f* D nlo 19 t x ) £* D x2 (o 2 .t 2 ) 
f+\ x {o x ,t x ) f* D 22 (o 2 ,t 2 ) 

f* D 91 lo lt t 1 ) f* D 92 lo 2 ,t 2 ) 



. . f* D 19 (o 9t t 9 ) 
f* D 29 (o 9 ,t 9 ) 

f* D 99 {o 9 ,t 9 ) 



3.12 



Each of the nine functions in each of the nine columns of 
equation 3.12 is defined in terms of the same independent 
variables, one of the nine stresses, and time. For example, the 
nine e n , 



3.13 



are defined in terms of the same independent variables, o y and t. 
These functions must also be coupled due to chemical bonding. 

When equations 3.13 are evaluated at some value of stress, 
say a, and at some time t, they degenerate to a set of constant 
values for strain, 

e D n (o lt t) - f* D n io 19 t) - z n 



The z n are the strains associated with the stress a, and time t. 
When evaluated at some other time, say t', these functions 
degenerate into another set of strain constants, 
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13 



where the z' n are the strains associated with the same static 
stress a, and time t'. 

One constant may be expressed in terms of another constant 
multiplied by another constant so the strain term z 21 may be 
expressed in terms of the strain term z 11 as, 



f* D 21 (o lt t) - z 2l - k 21 z lx - * 21 f*° 11 (a 1 , t) g 



and similarly for all other z n , i.e., z n = k n »f * D n (a 1# t) . The k n 



3.17 



jn are the scaling constants, and k n = 1. At time t' the strain 

!ii 

10 ig term z' 21 may be written as, 

ii 

M 

iu 

■ sn 

H and so on. 

Now, when t -> t' the values of f* D 21 (a 1# t) and f* D n (a 1# t) 
change, but the functions do not change. The scaling coefficients 
15 k' 21 and k 21 are independent of stress and time, so when t -> t' 
then f*° 21 (a 1# t) -> f* D 21 (a 1# t') and f*°„(a 1# t) -> f * D „ (a 1# t' ) . 
Equations 3.16 and 3.17 then yield, 

* 31 f^K, t') - Jr' ai f* D xl {a lt t f ) 

3.18 

K 21 K 21 
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for all time, and for all values of a r Following the same 
procedure for the remaining seven f*° n (a 1# t) will produce similar 
results and equation 3.12 reduces to, 

i f**ia<o a .t) ••• f* D 19 lo 9t t) 

k 91 f* D 92 (o 2 ,t) ... f* D 99 (a 9 ,t) 

where f* D u {a^,t) is a scalar valued function. Continuing this 
same exercise for each of the functions in the other eight 
columns in the matrix will then yield, 

1 i ... l 

^21 ^22 • " * ^29 

^91 ^2 " ' " ^99 

Because the nine constitutive functions f^jCa^t) are supposedly 
independent, each now represents a solution for a particular 
state of stress. Again, the simple proportionality between 
responses is an approximation and the result of imposing the 
small strain condition (see section 2) . 

All nine of the f*° %i ia } ,t) functions in equation 3.20 are 
scalar valued, and control all possible strains in all possible 
directions for any loading function operating in any given 
direction. However, equation 3.20 implies that each function is 
independent, but somehow associated with a particular type of 
stress, i.e., normal or shear, as well as a particular loading 
direction, where more that one function operates in a given 



3.19 



£*\ x (o x .t) 
f* D l2 (o 2 .t) 



f* D 19 (o 9 ,t) 



3.20 
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direction. Given the requirement for coupling in the strain 
components through bonding, and that loading from any direction 
produces deformation of the same set of bonds, this result makes 
no sense. Furthermore, the generalized non-linear functions of 
equation 3.20 must contain, as a special case, all possible 
linear viscoelastic solutions. So, equation 3.20 can not possibly 
be the final solution. 



3.4 The Final Reduction 

By using the reference frame transformation techniques 
already discussed the number of non-zero matrix coefficients 
associated with any degree of material symmetry or isotropy may 
be determined. 



Application of Normal Stresses 

Begin with the normal stresses, and set eight of the nine 
possible stresses equal to zero, with the exception of a,. 
Equation 3.20 then becomes, 



e D i (o 1 .t) - 



1 



f* D lx (a lt t) 



3.21 



k 91 



The shear strains are symmetric, so there are only six 
independent values for the k n and equation 3.21 reduces to, 
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e" 1 la 1 .t) - 



1 

k 61 
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From equation 3.11 the solution for total strain is, 



l 



*«1 
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Applying a 2 separately, and then applying a 3 separately produces 
similar solutions. Equation 3.20 then reduces to, 





1 


1 
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e D i (a jt t) - 
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f** 13 (o 3 , t) 
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where each stress is applied separately and alone , although the 
sequence of application is totally arbitrary. 



Multiple Planes of Material Symmetry 

Once again define two orthonormal coordinate systems, x f and 
z f , where the Xj*^ and Zj-Zj axes are coincident and define a 
plane, but where z 3 = - x 3 , the mirror image of the x 3 axis. 
Stress and strain in the x f reference system are a and e, while 
stress and strain associated with the z i reference system are a* 
and e*. Now consider the stresses and strains in the '2-3' and 
'3-1' directions. The direction cosines for this transformation 
are given by equation 2.27. Substituting those results into 
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equation 1.39 yields the vector transformation o* i} = o r} , and c*^- 
= € i} , except as follows, 



**33 " " °22* **31 - - <>31' e *23 - " 6 23 > € *31 " " 6 31 25 



None of the k^- coefficients associated with the normal strains 
5 in equation 3.24 can be eliminated, unlike the coefficients 
associated with the shear strains. From equation 3.20 the 
expression for the shear strain e 4 (..), when only a A is applied 
with respect to the xi reference frame, is, 

O 

\Q e 4 (a x/ t) - m 1A f* E (o 1 ) + * 14 f * D X1 (o 1# t) 

hp 3 - 26 

m 

-.7 : 

10 jg In the Zj frame the equivalent expression is, 

IB ' 

|3 6*4(0*!, t) - m 14 /**<o* 4 ) + *j 4 f** 11 (o* 1 ,t) 2? 

IU 

in 

P Flipping the z, axis and substituting the appropriate 

H 

transformations for stress and strain into equation 3.27 yields, 
-e 4 (o lf t) - jn 14 f*«(o l ) + ^ t f* D 11 (a l ,C) 2g 

15 Adding equations 3.26 and 3.28 produces, 

0 - 2in 14 f**(0 1 ) ♦ 2 * 14 /• D 11 (o 1 ,t) 3 2g 
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The constitutive functions are non-zero in value for any non-zero 
value of a t , and at any time t, so m u = k 14 = 0. The same result 
will be obtained for the transformation of e 5 , which then yields 
k i5 = *u = °' and equation 3.23 reduces to, 



€ i (a 1# t) 



l 

0 
0 



f* B {o x ) + 



^31 
0 

0 

*61 



f* D xx (o x ,t) 



3.30 



'3 Applying o z and a 3 separately, the same symmetry condition 
produces similar results and equation 3.24 becomes, 



= 1 3 



in 



in 



■"■21 

0 
0 



^32 
0 
0 



^33 
0 
0 



f* D xx (o lt t) 



3.31 



and likewise for the impulse response. 



10 Two/Three Planes of Symmetry 

Let's now discuss the case of two, and hence three planes of 
axial symmetry. In addition to the defined axial rotation of 
z 3 = - x 3 for a single axis of symmetry, the rotations for the 
condition of 2 axes of material symmetry now include z 1 = - x 1r 

15 z 3 = - x 3 , with z 2 = x 2 . The only non-zero direction cosines for 
these equivalences are given by equation 2.35. Equations 2.26 
yield the following transformations for £ 6 , and for the stresses 
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€ 12' °*1 " °1' 



o 2 , a* 3 - a. 



3.32 



Following the same above procedures will produce the equivalences 
k 6i = *62 = k 63 = °' and equation 3.31 reduces to, 



e D i (o j , t) 



1 


1 


1 


*21 


*22 


*23 


*31 


^32 


*33 


0 


0 


0 


0 


0 


0 


0 


0 


0 



3.33 



Normal Stresses and Complete Material Isotropy 

Begin with the results just obtained, where the two 
reference frames shall be initially coincident. Rotate the X| 
reference axes 90 degrees about the z 3 axis, such that the x t 
axis coincides with the z 2 axis, the z 2 axis with the -x, axis, 
and the x 3 axis coincides with the z 3 axis. Rotating the x t -x 2 
axes through some angle a about the x 3 axis produces the 
transformations from the x 1 -x 2 plane into the z,-z 2 plane given 
by, 

z 3 - x 3 , z x - XjCosa + x 2 sina, z 2 - -^sina + x 2 cosa 

3.34 

where a is equal to 90 degrees. From equation 1.39 the 90 degree 
rotation of axes yields the transformations, 
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°*1 " °2' °*2 " °1' °*3 " °3 

3.35 

e* x - e^, e* 2 - e 1# e* 3 - e 3 



Applying the stress a t in the x 1 axis direction the expression 
for strain in the x 1 direction yields, 



5 In the z 2 direction the expression for strain as a function of 
a* 2 is, 



O 
*0 



in 
in 

SB 



e D 2 *(a 2 *,t) - k 22 f* D 12 (o 2 *, t) 



rotate the z 2 axis so that it is coincident with the x, axis, and 
demand equivalence in the resulting strains (the condition of 

10 Si isotropy) . From the vector transformations equation 3.37 becomes, 

iU 



H <*V°2' t) - k 22 f* D 12 (o lt t) 

£ : 

Equating equation 3.36 with equation 3.38 yields, 

*22 f» B i 2 (°i.t) - '• D ii<«i» t) 

and the solution for the function f* D 12 (a 1# t) , 



3.38 



3.39 
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f* D 12 io lt t) - -A- f.^to^t) 



x 22 



3.40 



a, is a magnitude term, so the scalar function f* D 12 (a/t) scales 
in terms of the scalar function f* D n (<r # t). Multiplying through by 
the constant term l/k^ equation 3.33 becomes, 



e D i (o jt t) 



1 


l/*22 


1 


*« 


1 


*23 


*31 


^32/^22 


*33 


0 


0 


0 


0 


0 


0 


0 


0 


0 



/• t) 

f* D 13 {o 3 ,t) 



3.41 



Demanding the same degree of symmetry with regard to the Xj-Xj 
plane, equation 3.41 reduces further to, 



■^21 

*3i 
0 
0 
0 



l/*22 

1 

^32/^22 
0 
0 
0 



l/*33 

■^23/^33 

1 
0 
0 
0 



f* B 11 (a af t) 



3.42 



Redesignating the coefficients yields, 

f* D (a lr t) 





n X2 


«13 




"21 




^23 




«31 


^32 


"ll 




0 


0 


0 




0 


0 


0 




0 


0 


0 





f*°(o 2 ,t) 



f*°(o,,t) 



3.43 



where n t1 = 1, by choice of definition. 

Following the same procedures outlined in section 2 will 
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yield equivalence and symmetry in the non-diagonal coefficients, 
Skipping the algebra, equation 3.43 then reduces to, 





n 12 


n 12 




n n 




n !2 


n i2 




0 


0 


0 


0 


0 


0 


0 


0 


0 



f+ D (o 10 1) 



f* D (o 2 ,t) 



f* D (o 3 , 1) 



3,44 



where n 1t = 1. 



5 Shear Stress and a Single Plane/Axis of Material Symmetry 
^ Again set eight of the nine stresses to zero, except a t 

;3 Equation 3.11 then reduces to, 



iP e i (o 4 ,t) 

IS 
e 

%j where there are only six independent m i4 and k i4 due to symmetry 

iU . 

10 in in the shear strains. 

rf Let two orthogonal reference frames, the x f and z- systems, 

to be initially coincident. For a single axis of symmetry the 

axial transformations are again x t = z t , x 2 = z 2 , and -x 3 = z 3 . 

The solutions for the stress and strain transformations are 
15 a-j = a*|j, and c.jt..) = e* fj ( . • ) , except for those given by 

equation 3.25. From equation 3.45 the expression for the normal 

strain £,(••) in the x s reference frame is, 



^4 



f* B (o 4 ) + 



x 24 



3.45 
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e x (a 4 ,t) - f**(o 4 ) + ic 14 f* D 14 (o 4 , t) ^ g 



In the z f reference frame the equivalent expression is, 

e* 1 (o* 4 ,t) - m l4 f**(o* 4 ) + fc 14 f* D 1A lo* A . t) 47 

Substituting into equation 3.47 the appropriate transformations 
5 produces , 



13 
=0 



e 1 (o 4# t) - - th 14 f**(o 4 ) - * 14 f*° 14 (a 4 , t) 48 



Subtracting equation 3.46 from equation 3.48 yields, 

i'U 
■ •a 

W 0 - -2 ;n 14 f**(a 4 > - 2Jc 14 /• 1, 14 (a 4 , t) ^ ^ 

su 

Ifi Given that the constitutive functions are not equal to zero for 

o 

10 *2 any non-zero value of a 4 when time t > 0 then m u = k u = 0. 

Similar results will be obtained for the transformations of 
the other two normal strains and m 24 = = k 24 = k^ = 0. No 
information is obtained for the coefficients m a , m^, k u and k^, 
other than they are not equal to zero. These transformations 

15 however do yield = k^ = 0, and equation 3.45 reduces to, 
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e^o^ t) 



o 




o 


0 




0 


0 




0 




f**(o 4 ) + 










0 




0 



3.50 



Following the same procedures for the other shear stresses and 
resulting strains will produce similar results. 



Shear Stress and Two/Three Planes of Material Symmetry 

The axial transformations are again z 2 = x 2 , z 3 = -x 3 , and 
z 1 = -x r The associated direction cosines are given by equations 
1.40 - 1.41. This symmetry condition yields no new information 
pertaining to the coefficients m u or k^. However, for c 5 (..) in 
the x f frame we have, 



€5(04,0 - jn 54 f* B (o A ) + k SA f* D u(o A , t) 



3.51 



and in the z i frame, 

e* 5 (o* 4 ,t) - ji^ 4 f* E (o* A ) + Jc 54 f* D u (o* 4 , t) 



3.52 



Substituting the appropriate transformations for stress and 
strain into equation 3.52 yields, 



e s (o 4 ,t) - -/n^ f**(o 4 ) - * 54 f* c 14 (o 4 , t) 



3.53 



where f* E (-a) = -f* E (a) and f* D (-a,t) = -f*°(a,t), by definition, 
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Subtracting equation 3.51 from equation 3.53 produces, 



0 - -2 ^ f**(o,) - 2 k s< f*° 14 (a 4/ t) 



3.54 



Given that the constitutive functions are not equal to zero for 
any non-zero value of a 4 , then = = 0. Equation 3.50 then 
reduces to, 



0 




0 


0 




0 


0 

m u 


f**(o 4 ) + 


0 

K* 


0 




0 


0 




0 



3.55 



Continuing on with the other shear stresses, a 5 -> a 9 , will yield 
similar results for the other sets of coefficients. 

By rotating a reference frame 45 degrees any shear stress 
may be converted into a normal stress, although the magnitude 
differs. Given that f* D u (a,t) is scalar valued, the form of the 
function is unaltered by the reference frame transformation. A 
normal stress produces a normal strain, so f* D 14 (a,t) must be 
equivalent to f* D (a / t). Again, this should be obvious because 
simple rotations of spatial reference frame with respect to a 
body of material cannot cause changes in material properties. The 
same equivalence may be demonstrated for the remaining functions 
f*° 15 (a,t) -> f* D 19 (a,t). 



Complete Material Isotropy 

When the strain functions €,-(..) are evaluated for a 4 , or 
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for any other shear stress, no other non-zero values for the 
coefficients are found for the condition of complete material 
isotropy. However, we do find that m u = 1055 = m^, and that 

Given symmetry in the shear stress there is equivalence in 
the values of the constitutive functions when evaluated for each 
of these stresses, e.g., f* D (a ij# t) = f*°(a jl# t) . When the e f (..) 
are evaluated, six shear stresses yield no more information about 
the strain field than three shear stresses. So, equation 3.24 
becomes, 



e D i (o jt t) - 



1 


^12 


n i2 


0 


0 


0 


"12 


1 




0 


0 


0 


n 12 




1 


0 


0 


0 


0 


0 


0 




0 


0 


0 


0 


0 


0 . 




0 


0 


0 


0 


0 


0 







f* D (o lt t) 




f* D (a 2 ,t) 




f* D (o if t) 








f* D (o 5 ,t) 




f**(a 6 , t) 



3.56 



Inverting all the reference frame transformations and material 
symmetry or isotropy constraints will produce the original 9x9 
matrix of constant coefficients. However, there are only six 
equations and six unknowns, so equation 3.56 reduces to, 



e D i {a Jt t) 



"11 


"12 


"13 


"14 


" l5 


"16 


«21 


"22 


"23 


"24 


" 2S 


"26 


"31 


"32 


"33 


"34 


"35 


"36 


«41 


"42 


"43 


"44 


"45 


"46 


"si 


"52 


"S3 


"54 


"s5 


"56 


"61 


"62 


"63 


"64 


"65 


"66 



3.57 



where n u =1, by definition. The compete solution for strain 
becomes , 
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6,(0,., t) - lm u ] f* s (Oj) + [n^] f* D (o j9 t) ^ ^ 



As a general condition, symmetry in the matrix coefficients for 
any degree of material symmetry or anisotropy has not yet been 
demonstrated • 

5 3.5 Multiple Loads and the Non-Linear Solution 

The Stress Field and the State of Stress 

Following the discussion in section 2 define a Cartesian 
<5 reference frame, the x i system, where the nine a fj - act about some 
;D differential element of a solid material (Figure 2.2) . Shrink 

10 W that differential element to where SifSySz -> 0. The resultant 

W 

\U of the nine a % . stresses will then be the solution for the total 



state of stress at a point, a s (Figure 2.3). This resultant 



IS 

5 

M 

Hj vector has a magnitude and a direction of action with respect to 

!U 

jfi the X| axes which may be calculated. 

15 y h Now orient another orthogonal reference frame, the z i 

reference frame, with respect to a s such that one axis, say the 
z\ axis, is coincident and collinear with a s , and the z 2 " pZ 3 P lane 
is normal to a 8 (Figure 2.4). The z f system is now the frame of 
state for which the maximum state of stress is known, i.e., the 

20 type of stress, the magnitude, and the direction of action. 

Matrix Coefficients and Reference Frame Transformations 

Recall from section 2 that the number of non-zero matrix 
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coefficients is determined from various symmetry and isotropy 
conditions, and with reference to a primary Cartesian reference 
frame. The n ij|a and m ijkl matrices represent fourth order tensors , 
and like stress and strain tensors, the value of each coefficient 
varies according the orientation of respective reference frames. 
These coefficients transform from one reference frame, say the x- 
primary reference frame, to the frame of state, the z- system, 
according to the function 2.66. Again, the exception is complete 
isotropy, where the m* ijkl = m ijlcl , and the n* ijlcl = n ij|tl for all 
reference frame orientations. 

The Example of Isotropy and Non-Linear Viscoelasticitv 

Load an element of an isotropic, non-linear viscoelastic 
body such that the only stress a*, is oriented along the z, axis 
of the orthogonal z % reference frame (Figure 2.5a). Orient an x f 
reference frame initially coincident with the z f reference frame, 
then rotate the x f system about the x 3 axis (which remains 
coincident with z 3 ) so that z 1 -z 2 axes are oriented with respect 
to the x t -x 2 axes by some angle of rotation a (Figure 2.5b). In 
the x t -x 2 plane the vector a* t may be defined to be the resultant 
of two simultaneously applied normal stresses, a, and a 2 , 
operating coincident with the x 1 and x 2 axes. a* 1 now also 
represents the solution for the total state stress, a s , in both 
the z f and x s systems. 

For static loading the generalized governing constitutive 
equation of non-linear viscoelasticity is, 
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3.59 



which is valid for all orientations of orthogonal reference 
frames, regardless of the degree of material symmetry or 
isotropy. For complete isotropy the form of the matrix solutions 
are given by, 
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&12 


m x2 
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^2 


1 
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^2 


^2 
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0 




0 
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0 
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0 


0 


0 


0 


0 





3.60 



and likewise for the n^-, etc., for all possible reference frame 
orientations. 

From equation 3.59 the loading function a*, will produce 
three strains in the z i system (Figure 2.6a), €* 1 (a* 1/ t) = 
e* l1 (a* 1/ t), e* 2 (a* 1# t) = e* 21 (a* 1# t) and €* 3 <a* 1 ,t) = €* 31 <a* 1 ,t). 
The solutions for the strains e 1 (a* 1# t) and e 2 (a* 1# t) associated 
with the and x 2 axes are obtained by adding the appropriate 
vector components of t) and e* 21 (a* 1# t) (Figure 2.6b), 

e 1 (o* 1 ,t) - e* n (o* 1# t) cosa + e* 21 (a* x , t) sina 

3.61 

e 2 (o* 1# t) - e« u (a* a , t) sina + c* 21 (a* 1# t) cosa 

From equation 3.59 the solutions for the strains, €* t (a* 1# t) and 
e* 2 {a* 1# t) in the z- system are, 
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e* 1 (o* lf t) - e* xl la* lt t) - f* s (a* x ) + jj u /* i> (o« 1 , t) 

3.62 

€* 2 (o* 1# t) - e* 21 (o* 1# t) - £** {o* x ) + n 12 f* °(a* 1# t) 

Substituting equations 3.62 into equations 3.61 produces the 
solutions for €^(a* u t) and e 2 (a* 1/ t) in terms of a s = a* v 

e^o^t) - [m ll f* B (o* 1 ) + Ji n f* D <o* 1# t) ]coso + 
[ ^ f* B {a* 1 ) + n X2 f* D (o* r , t) ]sina 

3.63 

e 2 (a* 1# t) - [ % f^to^) + f+ D (Q* lt t) ] sina + 
[ m 12 f* B (a* 1 ) + n 12 f* D (o* ll t) Icosa 

So, as we saw in section 2, even with a simple isotropic 
material it is clear that only when f*(a,t) is linear, or may be 
approximated by a linear function, may the base reference frame 
solutions for strain be added in a linear fashion to obtain a 
correct solution. As with the non-linear elastic materials, a 
linear approximation of f*(a,t) will yield reasonable solutions 
for strain from the component stresses only if the solution for 
the state of stress does not exceed the range of accuracy of the 
linear approximation. 

Generalized Non-Linear Solutions for Static Loading 

To obtain explicit non-linear solutions for multiple loading 
conditions first define the stress components, the a- # with 
respect to a base Cartesian reference frame, call it x s . Next, 
solve for the magnitude and direction of action of state of 
stress vector with respect to the system. Now, orient another 
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Cartesian reference frame, call it the frame of state z f , such 
that one axis is coincident with state of stress vector, a s , and 
a 8 is normal to the plane defined by the other two orthogonal 
axes (Figure 2.7). Next, define the appropriate direction cosines 
5 for the system relative to the x. system, and compute the 

transformation solution for the two sets of matrix coefficients, 
the m^- and n^-, as given by equation 2.66. Finally, solve for 
strains in the z i system with the constitutive function, 

e*i<c,,t) - [a*^] f* B (o s ) + [n*^] f* D (o s ,t) 

3 • 64 

o 

10 ; B The index j may be 1, 2, or 3, and the state of stress vector is 
ID a normal stress vector aligned coincident with one of the z- 

y 

ffl axes. The index i will vary from 1 to 6 depending upon the degree 

W 

E of material symmetry or isotropy. 

Once the solutions for the €* i (<r s/ t) have been obtained in 

IU 

15 i.^ the z f frame of state they may then be transformed into solutions 
\ ' , for strain in x- base reference frame with the vector 
transformation equations. 

3.6 Time Limited Forms of the Non-Linear Solution 

In the z f frame the impulse response function f* E (a s ) may 
20 always be written as a power series expansion in a s , i.e., 
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f**(o a ) - a x o g + a 2 o s * + ... ♦ a n a s * - J] a n o," 

o-i 3 • 65 

X 

n-2 

where f* E (a s ) may always be expressed as the sum of a linear 
component, f* E L # <r s # and a non-linear component, f* E M L^ a s^ • 
The solution for the impulse strains in the frame of state, the 
€ *jt a s )' ma Y then be written as, 

3.66 

For linearity f* E (a s ) reduces to the first order term, with 
f* E L = f* E . 

At very large time, t->», the non-linear viscoelastic strain 
function e*(a s ,t->«) approaches an asymptotic limit for any 
constant a. At this time limit f*(a Q ,t->«) may likewise always 
be written as a power series expansion in c s (Figure 3.2), 

/j-i 3.67 

and expressed as the sum of linear component, f* L # ^ s , and a non- 
linear component, f* ML (^ 8 )- 

Now if the ultimate strain and the impulse strain may be 
expressed as power series expansions at t->«, then obviously so 
may be the delayed elasticity strain component, which is the 
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o 

£ i 

in 



difference between the ultimate strain and the impulse strain. 
The expression for total strain at t->* may then be written as, 



e*(o 8 , t-«) i - [m* n ]f* B (o 3 ) + ln* n ] £* D {o g , t-~) 



3.68 



or upon expansion, 



e* (a s , t-<*>) J - [m* n ] ^ a 3 + f*^(o s )] 



3.69 



The non- linear terms are the higher order expansions in a s . 
The solution for the linear strain term is, 



I 

b 

**j When the higher order non-linear terms are negligible then 

lOafi equation 3.69 degenerates to equation 3.70 and a 1 -> f* c and 

\t C % -> f* D (t->«) . 



3.7 Symmetry in the Coefficient Matrices 

The Work and Strain Energy Functions 

Figure 3.2 illustrates the one-dimensional form of the 
15 stress-strain relationship of equation 3.69. For either linear or 
non-linear viscoelastic materials the work, W, required to get 
from the undeformed state to the ultimate static deformed state 
at time t->« is given by the integral expression, 
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W - j o d € ( t-») 

0 



3.71 



This expression is the sum of an elastic potential energy 
function and a dissipation function, a heat term (the kinetic 
energy of deformation term, KE, is equal to zero) . When 
deformation is very slow and essentially isothermal then equation 
3.71 reverts to a strain energy function. 

A complementary work function, W*, may be defined by, 



W+ - j e ( t-») d o 



3.72 



and in a three dimensional format becomes, 



0 



J 

: , ( t-*co) da . 

3.73 



m - J ejU-co) dO] 



Rewriting equation 3.70 in a one dimensional format, ignoring 
heat loss, and substituting the result into equation 3.72 
produces a one dimensional recoverable strain energy function of 
the form, 

0 0 

W* - Je(t-«) do - j (e L -*-e liL )da 74 
oo 

a a 

- f e L da + | e^do - U* L + U*^ 

0 0 

where the total complimentary strain energy is the sum of linear 
and non-linear components. Expanding equation 3.74 into a three 
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!U 



dimensional format by substituting the generalized form of 
equation 3.70 into equation 3.73 produces, 



3.75 



In the z, frame of state there is only one stress function, a s , 
which is equivalent to a lf so equation 3.75 reduces to, 



w* - W* L + w*^ - f e L(i) do 8 + f e NLii) do g 



3.76 



P 
: n 

*5 The Linear Component of Work 

E 

if! Recall that in the z s frame of state the matrix coefficients 

in are transformations from the base reference frame, the 

i 

10= transformation being accomplished from equation 2.66. Obviously, 

M 

if the transformations are inverted we must obtain once again the 
same m ij|a and n ij|a . Furthermore, the products tt* ijkl *a 1 and n* ijkl »c 1 
|J of equation 3.70 will transform into m^^a, and n^^c, because a t 
and c, are scalar valued constants. An equivalent solution for 
15 the linear component of the strain field in the base reference 
frame may therefore be written in the form, 



where the o- } represent the stress components in the base 
reference frame. 



109 



The linear component of the strain energy function may now 
be defined in terms of the m iju# n ljia , a u and the component 
stresses, <7j. Because the strain field must be continuous and 
non-zero for any given applied stress field, the linear component 
of the complimentary strain energy function must exist. The 
following equivalence must therefore hold true, 

d*w+ L _ d*w+ L 
de i de j " de j de i 3.78 

and as before we find that the m^- = m }i , and n^- = n }i when j ^ i. 
So, the coefficient matrices are symmetric for all non-linear 
viscoelastic materials, for all degrees of material symmetry and 
isotropy. As with the linear materials the solution for symmetry 
is independent of strain, stress, the constitutive function, and 
system energy • 

3.8 The Non-Linear Constitutive Equations for Stress 

The constitutive functions for stress may also be expressed 
in terms of two matrices of tensor valued coefficients and a 
single scalar valued constitutive function, 

o(e jt t) i - [a ij ]f B (e j ) + [b ±j \ f D (e jt t) 

where the matrices of coefficients are symmetric in the indices 
ij, and a n = b n = 1, by choice of definition. The number of non- 
zero matrix coefficients for any degree of material symmetry and 
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isotropy are the same as for the m^- and n^- coefficients. Again , 
when written in the above format, the individual strains in 
equation 3.79 are constrained to be induced separately and alone. 

Now, we may define each component of strain in terms of a 
single vector, the strain state vector, e 8 , that defines the 
maximum magnitude of strain of a differential element of 
material. Following the above discussions we arrive at the 
function, 



10 ;D where the matrix coefficient are transformation determined 
tn through equation 2.66. The method of solving this equation for 



is 



the resulting stress field is the same as that used to solve for 
the strain field of equation 3.64. Once the a* f have been 

13 

h 4 determined in the reference frame of state, the o i may be 

IU 

15 m determined in the base reference frame from the transformation 



equations. 



3.9 Hereditary Integral Functions 

The creep compliance function f* D (a,t) and the relaxation 
modulus function f°(€,t) are always non-linear functions. They 
20 are non-linear functions in time, and may or may not be non- 
linear in stress or strain. Even when, for example, the function 
f* D (a,t) = f*(t)»a, the definition of a linear response, f* D (a,t) 
is still a non-linear function of time, although delayed elastic 
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strain, strain rate, and ultimate strain are proportional to a 
static o (Figure 3.3a). 

A non-linear delayed elasticity stress-strain curve defined 
by the function f* 0 (a,t=») may always be approximated by a series 
of piecewise continuous linear functions, as in Figure 2.1a. This 
implies that over the range in strain where each piecewise 
approximation is defined the functions f* D (a,t) may always be 
approximated in the linear form, i.e., f* D (t)»a. However, the 
various f*°(t) over the various linear approximation ranges need 
not necessarily be identical. Now, generalized function f* D (<r,t) 
defines the time dependent response for all possible materials, 
and includes as a subset of solutions the response of the linear 
materials. However, the linear approximations suggest that there 
is a possible subset of solutions where, for example, the delayed 
elastic strain, strain rate and the ultimate strain are 
proportional to some function of stress (Figure 3.3b), i.e., 

f* D (o,t) - f* D (t) *g(o) 3 ^ 81 

and similarly, 

t) - f D (t) - Me) 2 

These particular solution sets are suggested by the physics of 
deformation and the nature of the electromagnetic forces 
operating between atoms. 
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When the scalar functions g(a) and h(e) are approximated as 
linear functions over limited ranges in the stress-strain domain 
(Figure 2.1a) the linear hereditary integral functions may then 
be invoked as solutions to time dependent loading, aj(t), or time 
5 dependent deformation, 6j(t). This is convenient, for it avoids 
the complexity of calculating the stress and strain state 
functions, the orientation of the frame of state, the variation 
in the values of the matrix coefficients as the orientation of 
the frame of state varies with respect to the base reference 
10 frame, and the transformations of the resulting stress or strain 
vectors. However, the functions g(a) and h(e) are scalar valued, 
and scale the time dependent functions f* D (t) and f D (t) . So, we 

!H may defer from piecewise linear approximations and insert these 

W 

tfl functions directly into the linear hereditary integral functions 
IP 

15 * to obtain explicit solutions. Because there is now only a single 

S| loading function in the frame of state the resulting hereditary 

IU 

integrals for strain will take the form of, 

M 

€•,(■".<«.« - f**{o g {t)) + n*„ f g(o B { t') ) -2£L^LL dt f ^ ^ 

Again, the index j may be set to 1, 2, or 3. In general, the 
20 matrix coefficients m*^ and n*,-. also vary in time because the 

orientation of the frame of state, relative to the base reference 
frame, also varies. The one exception is, of course, the 
condition of complete isotropy where m i} and n f j are invariant for 
all orientations of reference frames. The solution for g(a) is 
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readily determined from Figure 3.1b. Obviously, an equation may 
be developed for stress which take a similar form. Developing 
that equation from the above discussions is straight-forward and 
warrants no further discussion. 

As mentioned, the deformation physics of some types of 
materials, such as the fibrous polymers, actually suggests 
equations 3.81 and 3.82 as a solutions for the time dependent 
components of the constitutive equation. Solid polymers, such as 
rubbers and plastics, are structurally defined by a chaotic 
lattice network of interlocking polymer fibers. As such these 
materials usually have relatively low densities because of the 
void space between fibers. Because of their structure polymers 
are often capable of very large strains (e.g., rubber bands). 
However, most natural crystalline materials, glasses, metals, 
etc., have a structural lattice framework defined by chemical 
bonds, and will only sustain strains on the order perhaps 1% 
before the bonds begin breaking. So, as a consequence of their 
structure we may anticipate that when the strains in the polymers 
exceed a few percent that changes in fiber geometries will 
produce non-linear changes in the internal distribution of 
stress, the components of which vary as trigonometric functions, 
and hence the constitutive functions. Obviously, in the case of 
large strains, the matrix coefficients need not be constant. 
Because large strains are probably more a function of fiber 
lattice deformation than bond stretching, the non-linear stress- 
strain relationships of the delayed elasticity functions are most 
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likely best defined by functions of the form given by equations 
3.81 and 3.82. The caveat being that strain magnitude be moderate 
enough that the stress-strain relationships between bonding and 
non-bonding forces still remain relatively linear, unless we can 
5 argue that equations 3.81 and 3.82 hold true for bond stretching 
as well. 

3*10 Temperature Effects 

The non-linear constitutive equations describing temperature 
induced stresses and strains were discussed in section 1 and 
10 *3 warrant no further discussion. Obviously, the constitutive 

; 0 properties of a material, will vary with changes in temperature, 
if! Following those discussions for the non-linear materials leads to 

y 

IB a temperature dependant solution of equation 3.58 of the form, 

m 

» 

O ei <e,o,,£) - [m^.] f* 8 (6,o.,.) + [n^] f* D ($.o jt t) 

\j 3,84 

!U 

in 
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SECTION 4 

ORTHOGONALITY AND UNCOUPLING 

Implicit to the derivations in the previous three sections 
are the assumptions that: 1) a volume of material under 
5 consideration was sufficiently large such that we considered a 
basic structural unit of a bonding framework between atoms or 
molecules; 2) that bonding geometries between atoms caused 
deformation in one direction to automatically produce deformation 
in the other directions, i.e., there is no orthogonality to 
10 ; -5 material structure or in the bonding arrangements, e.g., cubic, 
)0 orthorhombic, etc., or in the types of bonds - ionic, covalent, 

3 : 
s c= 

til etc; and 3) that deformation produced small strains, on the order 

W 

10 of a few percent at most, and negligible changes in internal 

IB 

lattice structure geometries and bonding angles. 

P 

ill 

15^ 4.1 Uncoupling and Orthogonality 

Orthogonality, Orthotropy and Normal Stresses and Strains 
The most obvious condition through which directional 
uncoupling in constitutive properties may occur is through 
bonding orthogonality, where the bonds between atoms are arranged 
20 at 90 degrees, and the loading is in the direction of the bonds. 
A few materials, such as table salt, sodium chloride (NaCl) , and 
lead sulfide (PbS) have orthogonal bonding (Figure 4.1). This 
bonding arrangement produces orthotropy, the condition of three 
axes of material symmetry. 
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Consider two simple three dimensional crystal lattice 
structures, cubic and orthorhombic, where all bonding is 
orthogonal (Figure 4.2) • Initially restrict the loading to be 
normal and coincident with the bonding directions, where each 
load is applied separately. Let the axes of an orthogonal 
reference frame be aligned and coincident with the bonds. First 
consider the cubic lattice, with identical bonds (Figures 4.1 and 
4.2a), and for the sake of argument let the constitutive 
responses 'perfectly' uncouple, eliminating the Poisson ratio 
terms, i.e., the lateral strains. Linear constitutive functions 
for the normal impulse strains would then be of the form, 



10 0 
0 10 
0 0 1 
0 0 0 
0 0 0 
0 0 0 



f * E • o 1 
f* B * a, 



4.1 



where the three f* E impulse functions are all equivalent and the 
Poisson ratio terms equal zero. There are no associated shear 
strains, c 4 -> e 6 , because of bonding orthogonality. Bonds need 
not be identical and the lattice structure may be orthorhombic 
(Figure 4.2b), so constitutive functions in each direction need 
not be equivalent. The general form of equation 4.1 is then, 
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Similar solutions may be obtained for the delayed elasticity 
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response and generalized linear viscoelastic functions for the 
normal strain responses become, 



e,(t) 



f*2' 
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1 



f* 2 °(t) • a 2 
f* 3 D (t) • Oj 



A solution for the stress functions would be of the form, 
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Equivalence may also occur between any or all of the impulse or 
delayed elasticity functions in equations 4.3 and 4.4. 

'Perfect' uncoupling probably never occurs, and for several 
reasons, the most obvious being that atoms in a crystal are not 
stationary, they always vibrate. However, a more compelling 
reason that 'perfect' uncoupling probably never occurs is that 
deformation of a bond entails the deformation of shared electron 
orbitals (Figure 4.3) and the associated electric and magnetic 
fields. The electromagnetic fields associated with the non- 
bonding/anti-bonding orbitals also provide reactive forces 
between the bonded atoms, and will also be effected by spacing 
changes between atoms (Figures 4.4 and 4.5). It is therefore 
unreasonable to expect or demand that distortion of the electron 
orbitals and their related electromagnetic fields be simply 
unidirectional. Likewise for those fields associated with the 
nuclei. Field distortion in the direction of the bonding will 
almost certainly result in lateral field distortions and produce 
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lateral strains. So, in general, a Poisson effect of some form 
must be admitted and equation 4.3 becomes, 
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Orthogonality. Orthotropy and Shear Stresses and Strains 

From previous discussion we know that the solution for the 
completely coupled linear impulse strain response for the 
condition of orthotropy is, 
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where a shear stress only produces a shear strain, and a normal 
stress only produces normal strains. This directional uncoupling 
in strains does not mean that the constitutive responses 
controlling shear deformation are uncoupled from those 
controlling the normal strains, it just means that they may be of 
different magnitude. 

Consider the cubic crystal lattice. Applying a shear stress 
to two opposite sides of one side of the cube (Figure 4.6a), the 
lattice racks, producing only shear distortions and no normal 
strains (obviously, only for very small shear strains) . Recall 
that for every a f j there is a corresponding equivalent a... The 
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constitutive function that describes shear distortion is, in 
general, not the same as that which defines normal deformation in 
the direction of the bonds (Figure 4.6b). When the lattice is 
deformed by an applied shear force, or equivalently by a normal 
force oriented 45 degrees to the bonding direction, more than a 
single reactive force acts to constrain deformation of the 
lattice. Both bonds, defined by the constitutive functions f t and 
f 2 , stretch simultaneously and define a coupled reaction. As bond 
angles change there is an additional reactive force, call it f 3 , 
established by the interaction of the electromagnetic fields of 
not only the bonding orbitals, but also of the non-bonding/ anti- 
bonding orbitals. Finally, there is a reactive force, f 4 , between 
the positively charged nuclei of the diagonal atoms that also 
acts to constrain lattice deformation. All these reactive forces 
act simultaneously. Simple stretching of the bonds in a normal 
direction does not bring into play all of these same force 
functions. 

Because there may be 3 different bonds there are 3 possible 
combinations of force functions, one for each of the 3 orthogonal 
planes defined by the crystal lattice faces (parallel faces are 
equivalent) . So, there are 3 possible shearing functions, each 
acting within one of the three orthogonal planes defined by the 
crystal lattice (Figure 4.2). The generalized impulse response 
then becomes, 
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while the delayed elasticity strain response is defined by, 
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The three normal functions may be equivalent, and the three shear 
functions may be equivalent. However, the shear and normal 
functions are most probably never equivalent, and for the reasons 
just discussed. 

Figure 4.7 illustrates an orthorhombic crystal lattice 
oriented so that one axis, the x' 2 axis of the x' 5 reference 
frame, runs diagonally through two corners of the structure. 
Applying a load coincident with that axis causes the entire 
lattice to deform, stretching all bonds and changing all angles 
between the bonds. All forces acting between the atoms come into 
play simultaneously and produce a coupled response. The 
constitutive response in this particular reference frame is 
therefore defined by a single scalar valued impulse function and 
delayed elasticity function. 

The static load, a' 2 , may also be defined as the resultant 
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of as many as six separate loads, the a jf operating in the 
reference system defined by the direction of bonding, the x f 
reference frame, where these loads now operate simultaneously in 
the direction of the orthogonal bonding and/or in the orthogonal 
planes defined by the lattice. It is evident that unless loading 
is normal to, and coincident with, the direction of bonding, or 
shear loading is directionally constrained to be within one of 
the planes defined by any two orthogonal set of bonds, the 
constitutive response will be defined by a single scalar valued 
function. Partial uncoupling only occurs when loading is 
constrained to be coincident with the direction of bonding, or 
within a plane defined by two orthogonal sets of bonds. Hence, 
orthogonal bonds and lattice planes define directional 
singularities in the constitutive response. 

4*2 Orthotropy and Symmetry in the Matrix Coefficients 

The Linear Impulse Response 

Recall that the test for symmetry in the matrix coefficients 
involves the use of a strain energy function, U, and the 
equivalences , 

dU _ dU dU _ dU 

da j do j do j do i de i de^ de^de^ 4.9 

Applying this test to the impulse response of equation 4.7, for 
the stresses a t and a 2 , produces the equivalence, 
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^12 f*2* * ^1 *V 



4.10 



Given that f*^ ^ f *2 E ' then m 12 ^ m^, and so on. So, the matrix 
coefficients of equation 4.7 are not symmetric, at least in the 
format in which the equation is written. 

However, one constant may always be written in terms of 
another constant so equation 4.7 may be rewritten as, 
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where the f*j E = d^f*, 8 . Multiplying through by the dj's yields, 
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where m' u = d-mi;-, and d, = 1, by choice of definition. Equation 
4.12 now has the appearance of a completely coupled set of linear 
elastic functions. 

Substituting the equivalences of equation 4.12 into equation 
4.10 (where the appropriate number of a- are applied) produces, 
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where m' 12 = m' 2 i- Likewise for the other m'... So, regardless of 
whether the constitutive functions are completely coupled, or 
partially uncoupled, the linear impulse response may always be 
written in terms of a symmetric coefficient matrix and a single 
scalar valued constitutive function. 

The Delayed Linear Elastic Response 

The six creep compliance functions of equation 4.8 are non- 
linear functions of time. However, when evaluated at t->« the six 
f*j D (t->«>) are simple constants. Letting the f*j D <t->«>) = 2. the 
Zj's may be expressed as z } = c^z,, where c 1 = 1. Equation 4.8 now 
becomes , 




4.14 



which is of the same form as equation 4.11. Multiplying through 
by the Cj's produces a set of matrix coefficients, the n'f-'s, and 
a solution of the same form as equation 4.12. Invoking the same 
test as applied to equation 4.12 will show that the n'.j 
coefficient matrix is likewise symmetric. 

It is now a simple matter to rewrite equation 4.8 for all 
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# # 



time t. Factor out the value of unity from each f*j D (t) function 
of equation 4.8, where unity for each function is defined as 
Cj/Cj and each time independent Cj is as determined from equation 
4.14. Multiplying each column of the matrix by the appropriate Cj 
produces matrix symmetry and equation 4.8 becomes, 
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where ff'^Ct) = f* j (t)/c j . At t->» the f'*j<t->«) -> z,. 

Partial Orthogonality and Uncoupling 

There are three different sets of bonds in the lattice 
system of Figure 4.8. Each may be a different type and/or between 
atoms of different elements. Only the bond sets between the atoms 
labeled A and B, and C and D are directionally orthogonal. 
Because of potential differences in bonding, it is necessary to 
allow the three sets of bonds to be defined by different 
viscoelastic functions, call them f,, f 2 and f 3 . In addition to 
these forces there are those reactive forces that arise between 
the electromagnetic fields of the bonding and anti-bonding/non- 
bonding orbitals, f 4 and f 5 . Finally, there are the reactive 
forces between the nuclei of the atoms B and C, the force f 6 . 
Considered separately each of these functions has a potentially 
different time dependent response (see section 5) • 
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Now apply a constant static stress a 1 in the x, direction. 
This load propagates through bond sets 1 and 2 into bond set 3, 
producing the time independent impulse response. Is it reasonable 
however, to presume that time dependent deformation of some 
5 combination of bonds, say the B-C and C-D bonds, may possibly 

halt prior to the point in time when deformation of the A-B bond 
halts? 

When the bonds between the D and C atoms are stretched there 
is a corresponding change in the orientation and intensity of the 
10 electric and magnetic fields operating between the atoms, defined 
-3 by both the bonding and anti-bonding orbitals and the nuclei. 
This reaction produces the coupled constitutive responses f 2 and 

tn f,. These changes in field orientation and distribution produce a 

W 

19 displacement in atom C, which produce a simultaneous change in 

5 T"S 

15 s. the fields between atoms C and B, establishing the reaction f 3 . 

n 

*iJ This leads to a displacement in atom B, which establishes the 

Ifl reactions f 5 , f 6 . Any change in the spacing between atoms B and C 

p 

i2 causes a corresponding reaction between atoms C and D and the 
constitutive reactions f 2 and f 4 . Simultaneously, the change in 

20 the spacing between atoms B and C leads to a displacement in atom 
B, which establishes the reaction f v and displacement of atom A, 
and so on. So, it makes no difference whether the applied force 
a t produces time independent or time dependent displacements, the 
reactions are coupled in both space and in time, and the Poisson 

25 ratio terms are time independent coefficients. 
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Other Examples of Orthotropy 

Three dimensional bonding orthogonality implies material 
orthotropy, but orthotropy does not imply orthogonality. One 
example of material orthotropy without bonding orthogonality is 
the face-centered cubic bonding system, common to many metals 
(Figure 4.9). This type of bonding arrangement produces a lattice 
with three axes of identical symmetry. In the pure metals, 
composed of a single element, the bonds are all equal and 
oriented at 60 degrees. From this example we may conclude that 
the constitutive properties of materials characterized by non- 
orthogonal bonding arrangements will be defined in terms of a 
single scalar valued constitutive function, whether loading is 
coincident with the axis of material symmetry or not. 

4.3 Uncoupling and Other Material Symmetry Conditions 

Orthotropy is just one of three material symmetry conditions 
where the shear stresses produce only shear strains, and normal 
stresses produce only normal strains. The other two symmetry 
conditions that display such uncoupling are: one plane of 
isotropy and one axis of symmetry; and complete isotropy. The 
constitutive properties of isotropic materials have already been 
discussed in detail and partial directional uncoupling is not 
possible for this symmetry condition. 

It is not necessary that two or three sets of bonds be 
orthogonal for partial uncoupling in the constitutive response to 
occur. One set of bonds orthogonal to a sheet of atoms that have 
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no orthogonality may also produce uncoupling and material 
orthotropy. Furthermore, it is possible that a given material may 
be compositionally and/or structurally laminated in one or more 
directions, where each layer has a distinct and separate chemical 
5 composition and/or structure. Each such layer may be thick enough 
and structurally complex enough to have its own distinct 
constitutive behavior. 

One Plane of Isotropy and One Axis of Symmetry 

Like orthotropy, this symmetry condition does not 
10 i3 necessarily imply orthogonality. However, orthogonality is 

=P possible. There are few naturally occurring minerals that have 

If! laminated structures, one of which is graphite (Figure 4.10). 

W 

fn There are even fewer lattice structures that can produce 

/ crystalline isotropy. One such structural arrangement is the 

15 l*j hexagonal benzene ring, a sequence of carbon atoms that share 

ilj 

jp resonating hybrid double-single covalent carbon-carbon bonds 
H (Figure 4.11). Within the plane of the benzene ring there are 
three axes of associated symmetry, oriented 60 degrees to each 
other. This type of symmetry produces material isotropy in the 

20 plane of bonding. The carbon atoms in graphite are arranged in 
sheets of interconnected benzene rings bonded together by weak 
van der Waals bonds oriented 90 degrees to the sheets of carbon 
atoms. Given the orthogonality of the Van der Waals bonds to the 
carbon sheets it is clear that there will be directional 

25 differences in material constitutive properties between the axial 
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direction coincident with the van der Waals bonds and those in 
the plane of the carbons sheets. 

Let's now define a reference frame for graphite such that 
the x 1 and x 2 axes are oriented within a plane defined by the 
carbon sheets, and the x 3 axis is coincident with the van der 
Waals bonds. The constitutive equation that defines the 
singularities in the delayed elasticity response for this frame 
of reference are then, 
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Because of isotropy within the x 1 -x 2 plane, there are only three 
different delayed elasticity and impulse functions. The function 
£'**(t) defines the isotropic properties within the plane of a 
carbon sheets, f'* 2 D (t) defines the constitutive response normal 
to the carbon sheets, and f'* 3 D (t) defines the shear response 
function in the two planes normal to the plane of the carbon 
sheets. Using the same scheme as applied to the orthotropic 
material we may show that the matrix coefficients are symmetric. 

Complete Anisotropy 

Aside from the condition of complete isotropy, those 
symmetry conditions that do not produce uncoupling of the normal 
and shear responses are complete anisotropy, and the condition of 
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one axis of symmetry. No matter how a reference frame is oriented 
with respect to a body of an anisotropic material there are no 
axes or planes of symmetry. This complete lack of symmetry 
produces the impulse response function, 
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The matrix coefficients are all non-zero in value, so any applied 
normal stress produces not only three normal strains, but also 
three shear strains in the three axial planes. Likewise any 
applied shear stress produces three normal strains and shear 
strains in all three axial planes. Any applied stress causes not 
only normal deformations, but also racking of the lattice 
structure in all three directions. Hence, partial directional 
uncoupling in constitutive properties is not possible. 

One Axis of Material Symmetry 

When a material has one axis of material symmetry the 
impulse strain response will be given by the function, 
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Any normal stress produces not only normal strains, but also a 
shear strain in the 1-2 plane, i.e., the strain e 6 . When a shear 
stress is applied the shear strains do not completely uncouple. 
However, the deviatoric racking of the lattice structure is 
5 directionally dependent upon the applied load, and not all shear 
stresses produce normal strains, e.g., a 4 and a 5 . So, partial 
directional uncoupling in the constitutive functions is possible 
for this symmetry condition. 



Hydrostatic Loading and Singularities Due to Structural 
10 f 5 Symmetries 

t.ri 

^0 Consider once again the cubic crystalline structure of 

H 

in common table salt (Figure 4.1). For a crystal of this material 

y 

i'O there will be no racking of the lattice structure under 
hydrostatic loading. All net resultant forces act in the 

15 ^ direction of the bonds. Hence, hydrostatic loading also produces 
a singularity in the constitutive response. However, because the 
\i structure of the crystal is cubic and all bonds are identical and 
normal to one another there is no difference in the constitutive 
response as compared to equations 4.5, 4.7 and 4.8. 

20 Now consider the body centered cubic crystalline structure 

of a material like pure metallic iron (Figure 4.12). Clearly, 
there is no orthogonality to bonding for this type of crystalline 
structure. However, all bonds are identical and the bonding about 
the central atom is geometrically symmetric. So, hydrostatic 

25 loading cannot cause racking of the crystalline lattice structure 
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and the resulting deformation is the same as if loading had been 
applied only in the direction of the bonds. Hence, for 
hydrostatic loading this type of crystalline structure also 
produces a singularity in the constitutive response of the 
material. 



4.4 Non-Linear Materials and Matrix Symmetry 

Matrix Symmetry and the Non-Linear Impulse Response 

Again assume orthogonality in bonding, resulting in material 
orthotropy. Let the constitutive response to loading be non- 
linear for even very small strains. The most general solution for 
the resulting impulse response functions is then, 
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where the individual stresses are applied separately and alone. 
Each non-linear constitutive function may be expressed in terms 
of a power series expansion, and each expansion may be expressed 
as the sum of a linear and nonlinear function, e.g., 



f*i B (°i) ■ a i °i + t>i a i 2 + • • ■ 
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NL 



and so on. Rewriting equation 4.19 in terms of the non-linear and 
linear responses produces, 
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The linear term of equation 4.21 is identical in form to 
equation 4.7. So, rewriting the aj in terms of a 1 and a series of 
scaling constants, Cj, produces, 
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where aj = cfa^. Multiplying through by the Cj's yields, 
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where in this format the m'^. coefficients are symmetric. 
Now, factor out unity from each of the six non-linear f*j E (aj) 
functions, where unity for each function is defined as Cj/Cj, as 
determined from equation 4.22. Equation 4.19 then becomes, 



c i/ c i • f*i*(<*i) 

c 2 /c 2 * £+ 2 B (o 2 ) 

C3/C3 • f* 3 *(a 3 ) 

c 4 /c 4 • f* 4 *(a 4 ) 

c s /c 5 • f* 5 *(a 5 ) 

cJ c 6 * f*6*(°6> 



1 


^2 




0 


0 


0 




1 


**23 


0 


0 


0 


^1 


^2 


1 


0 


0 


0 


0 


0 


0 


1 


0 


0 


0 


0 


0 


0 


1 


0 


0 


0 


0 


0 


0 


1 



4.24 



133 



Multiplying the by the appropriate Cj, and defining 
f'*j E (<Jj) = f*j E (<7j)/Cj produces the symmetric matrix solution, 
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Matrix Symmetry and the Non-linear Creep Compliance Response 

Given the condition of orthogonality and resulting material 
orthotropy, the most general non-linear form of a delayed 
elasticity function is similarly, 
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At large time, t->», transient time dependent deformation ceases 
and all f * j D (t->« / a J .) become constant valued, independent of time, 
and may be expressed solely as functions of stress. The six 
f *j D (t->« / a J -) may then be rewritten as power series expansions of 
stress, and equation 4.26 takes the same form as equation 4.19. 
The same procedures applied to equation 4.19 may be used to show 
that the n-. matrix may always be expressed as a symmetric 
matrix, producing a solution for equation 4.26 of the form, 
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where f'^t,^) = f * j D (t / a j ) /d j# and the n' fj = n^d,-. d, = 1, by 
definition. 

These same procedures may be applied to the other conditions 
5 of material symmetry / as with the linear materials, so no further 
discussion is warranted. 

J8 4.5 Induced Singularities 

ffi it is also possible to induce directionality in manufactured 

(0 materials. Reinforced fiber composite materials, for example, are 

10 = constructed by imbedding fibers or a woven cloth of one type of 

(3 

material in a matrix of another type of material. Typical 

IU 

jfs examples of such types of materials are fiberglass and carbon 

d fiber reinforced plastics, steel belted rubber tires, steel 

i,™ 

reinforced stressed and unstressed concrete, etc. The purpose of 
15 imbedding fibers in a matrix of another material is usually to 
limit tensile deformation and provide a degree of tensile 
strength that the matrix alone does not possess. The coupling 
mechanism between the imbedded fibers and the matrix is most 
often mechanical, but may be chemical as well. 
20 Now non-reinforced matrix materials are, like rubbers and 

plastics and concrete, typically isotropic. If the imbedded 
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fibers are direct ionally oriented, as in a 90 degree blanket 
weave, they will impart a directionality, or even a directional 
singularity, to the constitutive response of the composite 
material. A 90 degree blanket weave will produce material 
5 orthotropy (Figure 4.13), with identical properties within the 
plane of the weave, provided the type and number of fibers in 
each direction is identical. A directional singularity normal to 
the blanket weave is also possible, depending upon how the woven 
fiber blankets are laminated, and so on for the planes of shear 
10 normal to the axes of singularity. Using a different number of 
fibers in one direction compared to another, or imbedding 

■? 

parallel fibers in just one direction produces not only 
j.fj orthotropy, but also a directional singularity in the axial 

!S direction parallel to the fibers, and so on. Using fibers of 

fQ 

15 i different material properties, like carbon and steel, oriented at 

O 

90 degrees, would likewise produce directional singularities in 

!U 

in the axial directions parallel to each type of fiber. 

w 

Constitutive Property Approximations 

It is not uncommon for the constitutive properties of the 

20 reinforcing fiber material to be markedly different than that of 
the surrounding matrix material. All the constituent materials 
are always fundamentally viscoelastic. However, the disparity in 
the delayed elasticity response of one material in comparison to 
the other may be such that the time dependent response for one 

25 may dominate to such an extent that the time dependent response 
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of the other may be ignored in defining a constitutive function 
for the composite material. Steel belted rubber tires and 
reinforced concrete would be examples of such possible types of 
composite materials. The delayed elasticity responses for rubber 
and concrete are far more pronounced than they are for steel (the 
delayed elasticity response of steel is insignificant enough that 
most engineers are unaware of the existence of this component 
response, so it is ignored for virtually all applications) . 

Assuming linearity in the constitutive responses for both 
the fibers and matrix materials the governing constitutive 
response for the fibers may be written as, 

e(t) iif) - lm ijif) )f* B {f) + [n ijif) ]f* D (t) {f) 



and for the matrix, 

e<t>i(») - + l**m*lf* D MM 4 29 

If the creep compliance response of the fibers is 
negligible, then regardless of the relative volume of the fibers 
with respect to the matrix material we may be able to ignore that 
response component, and the resulting constitutive function for 
the composite material would be of the form, 

€(t) i(c , - [m ijic) ]f** ic) + ln ijic> ]f* D (t) {a) 
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If the matrix material is isotropic then a 60 degree weave 
would produce isotropy for the composite material in the plane of 
the weave, while a 90 degree weave produces orthotropy. Needless 
to say, a similar scheme may be applied to non-linear materials 
as well. 



4 • 6 Temperature 

The constitutive response of any material will undoubtedly 
vary with temperature. Letting 6 represent temperature, equation 
4.15 may be rewritten as, 
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Likewise for the linear impulse response, and for the non-linear 
materials. The only demand applicable to these functions is that 
there be no temperature variation within the volume of material 
under consideration. So long as strains are small and changes in 
lattice structure geometry are small and negligible, the matrix 
coefficients should show no change in value. 



Temperature Induced Strains 

If a volume of some mechanically unstrained and unstressed 
material changes temperature it should either contract or expand. 
For an unstressed volume of material this produces a constitutive 
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function of the form, 



4.32 
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where 60 = 0 - 0 is different than the initial reference 
temperature 0 O . The strain functions A;(t,£0=O) = 0, by 
definition. The time variable defines the time required to reach 
strain equilibrium from the point in time of the chanqe in 
material temperature from the reference temperature, 0 Q . 

Thermal expansion and contraction will stretch and compress 
bonds between atoms, like an applied mechanical stress. So, the 
A ; (t,£0) must, in general, be characterized similarly, i.e., 



The A E ,(£0) and the A D j(t,$0) are coupled functions, so for 
completely anisotropic materials equation 4.32 reduces to, 
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15 and so on for the other degrees of material symmetry. However, as 
with mechanical strain functions, bonding orthogonality will 
produce separate thermal strain responses. Consider, for example, 
the condition of material orthotropy, where bonding is orthogonal 
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in three axial directions. Equation 4.32 now reduces to, 
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All three elastic and transient functions may be different, or 
may be equivalent, or any combination of any two functions may be 
5 equivalent, depending upon the degree of orthogonality and 
symmetry . 
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SECTION 5 

FLUIDS, UNCONSTRAINED DEFORMATION* AND INELASTIC DEFORMATION 

Unlike the solids the fluids are materials where there is no 
chemical bonding between molecules, or groups of molecules. 
Deformation is therefore potentially unconstrained for any 
magnitude of applied deviatoric stress. 

5.1 Response of 'Simple' Fluids to Deviatoric Loading 

A 'simple fluid' shall be defined as a material that can not 
support an applied shear stress of any magnitude, where the lack 



13 

10 J y of bonding or any other forces of significance acting as coupling 



^ agents between molecules, or clumps of molecules, precludes any 

lis 

|.9 resistance to shear forces of any magnitude, other than a 



frictional resistance. The relationship between a shear stress, 



h '4 <* u , and deviatoric (shear) strain rate, 3€' ij /3t, or deviatoric 

fU 

15ff] rate of deformation, D'^, may then be defined as, 



or 



where rate of deformation, D {J -, is defined as, 
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The V| are velocity terms, and the x f are spatial coordinates. 
The rate of deformation function, D ij# like shear strain, is 
always symmetric, and D'^- = D'j,-/ when i ^ j. All sorts of linear 
5 and non-linear relationships are possible for these functions 
(Figure 5.1). As we shall see from the following discussion, 
equations 5.1 and 5.2 are not even approximately equivalent 
unless total deformation and total strain are small, on the order 
¥ of a percent or so. Unlike the strain rate function the rate of 
10 lB deformation function is always linear, even for very large 

-M deformations. 

W 
IB 
ip 

5.2 Small Strain Rates and Rate of Deformation 

O 

From strain theory we know that we may write the expression 



k -4 



iU 

if\ for finite strains as, 



u 



_1 dux duj^ i_ ( du k 

15 6 ^ 2 1 dxj + 3x/ + 2 { dxj dx/ 5.4 



where the u f , u j# and u k are displacements and the x k and x. are 
the material coordinates. For the small strain condition the 
higher order partial derivative terms are much smaller than the 
first order terms, so equation 5.4 reduces to, 
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Now, the x s , etc., of equation 5.5 may be also taken to be 
spatial coordinates for small strains. So, differentiating 
equation 5.4 with respect to time produces, 





5.6 



2 dt dx. 



2 dx j 




where du^dt = v s and the v { are velocity terms, etc. Hence 
equation 5.6 becomes, 



which is equivalent to equation 5.3. So, the small strain 
restriction on all possible forms of equation 5.1 for approximate 
equivalence with equation 5.2 should now be clear. 



The simplest of all possible 'simple' fluids are Newtonian 
fluids. By definition, Newtonian fluids have a linear 
relationship between stress and the rate of deformation. The 
proportionality constant between stress and rate of deformation 
is the 'viscosity' of the fluid, a pressure and temperature 
dependent measure of resistance to deformation under a shear 




5.7 



5.3 



The Definition of Newtonian Fluids 
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stress. All gases behave as Newtonian fluids, do many liquids. 
From equations 5.1 and 5.2 the constitutive equations of 
Newtonian fluids may be defined as, 

°* - ^St* 5.8 



or 



By definition, M is dynamic 'viscosity'. The factor '2' is a 
constant of convenience. These algorithms are valid only for 
dynamic flow conditions. They do not define the static or dynamic 
behavior under hydrostatic loading. 

5.4 The Continuity Condition of Mass Flow 

The Mass Balance Equation 

Consider a reference volume (Figure 5.2) located with 
respect to some spatial reference frame x-y-z, and through which 
mass may flow in all directions. The rate of mass, p, accumulated 
per unit of time, dp/dt, in the reference volume Ax*Ay*Az, is the 
difference between the mass flowing into the volume and the mass 
flowing out. In the x direction, the mass flowing into the 
reference volume per unit of time, at velocity v x , through the 
area Ay*Az at reference point x, y, z is given by Ay*Az*p # v x . The 
mass flowing out of the reference volume per unit of time, and 
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through the area Ay*Az located at reference point x + Ax, y, z is 
given by Ay # Az»[p*v x + A(p*v x )]. The difference between these two 
expressions defines the mass rate of change, 

Ax Ay Az -fl - Ay Az p v x - Ay Az[p v x + A (pv x ) ] 

5.10 

Similar expressions may also be derived for flow into the 
reference volume from the y and z directions. Summing these 
expressions and taking the limit as Ax, Ay, and Az -> 0 produces 
the following differential expression for the time rate of change 
in mass density, 

3£ r d{pv x ) d{pv y ) d{pv z ) 

dt " dx dy dz 5.11 



Reference Frame Transformation 

Performing the indicated differentiation on the right side 
of equation 5.11 produces, 




5.12 



The left side of this equation is now the time derivative of 
density along the path of fluid motion, the 'substantiative' or 
material derivative. The right side is the solution for the local 
change in density about some moving reference molecule as a 
function of the divergence in the fluid velocity field, V«v, 
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where by definition, 



_ dv r dv v dv. 

V • v - - + -i- + £ 

3x dy dz 5.13 



Equation 5.12 may be rewritten as, 



4£ - - p (v * v) 

5 . 14 



where , 



i?t at x a* r a* 5.15 



Equation 5.14 is the continuity equation of compressible 



fB fluid flow. If a fluid is nearly incompressible, or when density 

^ changes during flow are negligible, then Dp/Dt = 0 or Dp/Dt » 0 

S3 

10 N and, 

!U 



V • v - 0 

5.16 



For static conditions the v | = 0, and equation 5.15 reduces to 
the obvious equivalence. 

5. 5 Density, Volumetric Strain and Rates of Deformation 

15 Let's now develop some relationships between material 

density, strain, strain rates and rates of deformation that we 
shall use in upcoming discussions. 
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Let an initial reference volume, as in Figure 5.2, have 
initial dimensions and volume of unity. Let the finite strains, 

represent the changes in the length of each side of this 
unit volume. Let V c represent the initial unit volume, V f the 
strained volume, and e v the volumetric strain. The solution for 
the volumetric strain may now be written as, 

(l + e ir ) (l + e 22 ) (l + e 33 ) - V Q 
6v " v 0 5.17 

Given that V Q = 1, then the time derivative of e y is, 
d It - ~4z I (1 + e xl ) (1 + € 22 ) (1 + e 33 ) - 1 ] 



dt dt ^ 1X ' * ' ' 5.18 

For very small strains the higher order terms of equations 5.17 
and 5.18 may be ignored yielding the solution, 



d e v m d ( x ***** m D 



where, by definition, 



D • 2 U 



Equations 5.14 and 5.19 then yield the following equivalences and 
approximations , 
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JcJr 



|J - - p (V-v) - - p^ - - p dfc r dt 5>21 

Now, the density, p, of any material at some initial 
volume, V 0 , and another volume, V 2 , and density condition, p 2 , may 
be written as> 

m m 
9 ' ^ ' P2 ^ 5.22 



where m is mass, a constant, and where V 2 = V 0 + dV. e v represents 
the differential volumetric change of the initial unit volume of 



U material V Q . So, from the definition of V 2 it is obvious that we 

in 



may write, 



v 2 - v 0 dv 



1<L v o 1 " dV " * v 5.23 

u 
IU 

)H Defining p 2 in terms of the differential change in density, dp, 
gives us, 



dp 

5.24 



From the definitions of p 2 and V 2 the second of equations 5.22 
15 may be rewritten as, 

+ " ( V.?*,) 5.25 
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Substituting the solution for m from the first of equations 5.22 
then yields, 



p v Q 

p + dp ( v Q + dv ) 5.26 



or upon rearrangement where V = 1, 



Wrt A dV 

dp - -p 



< 1 + dv) 5,27 

which holds true for any material. Confining interest to small 
volumetric changes, where dV « 1, equation 5.27 reduces to, 

dp " " p (ifdV) ' " P dV • " p m " p 6 " 5.28 



5.6 Hydrostatic Loading and the Kinetic Equations of State 

No known fluid demonstrates unconstrained volumetric 
deformation under a constant hydrostatic load. Unlike the 
deviatoric response, the degree of deformation associated with 
any given hydrostatic stress is constrained, regardless of the 
magnitude of the applied stress. As with the solids, it must be 
that the electromagnetic and thermodynamic forces acting at the 
atomic level prevent the collapse of molecules in fluids onto 
each other. These two markedly different responses in the simple 
fluids demand that each response be described by separate and 
uncoupled constitutive functions, where uncoupling is a 
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consequence of the lack of bonding or other coupling forces 
acting between molecules. 



Definition of the Stress Functions 

Experimental observation tells us that for the simple static 
no-flow condition, or the uniform flow condition, a simple fluid 
will not support any applied shear stress. The total state of 
internal stress is then purely hydrostatic. Total hydrostatic 
pressure P is defined as, 

r„ - -f«. - - Tkk - - Tl1 + Tzz + r " 



— 3 5.29 

where T lM is total stress, and P is the mean hydrostatic 
pressure. The (-) indicates compression. 

Consider steady state flow, where a velocity gradient is 
permitted within a flow stream. Let some simple fluid, a liquid, 
be flowing through some channel (Figure 5.3a). Viscosity and 
frictional drag at the fluid-solid interface produce the velocity 
profile internal to the flow stream. Let flow rate be constant 
and the internal velocity profile be constant, with gravity as 
the driving force. Similarly, let gravity alone produce the 
hydrostatic stress. Under these conditions the hydrostatic stress 
remains unchanged in the flow column because there is no change 
in fluid density or volumetric flow rate. This condition demands 
no divergence in velocity field, so from equation 5.21 D,^ = V*v 
= de v /dt = 0. The hydrostatic stress, which is in equilibrium, 
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must vary with column height, where T fi = -P = -p*g*h. h is fluid 
column height, and g is gravitational acceleration. Again, the 
{-) simply indicates compression. 

Total stress, P, may be expressed in terms of two separate 
components, an applied hydrostatic stress, a H , and a 
thermodynamic stress (pressure) p, a function of the 
thermodynamic equation of state. At equilibrium the total 
hydrostatic stress P equals the thermodynamic stress p. 

The thermodynamic stress p may be defined by some material 
dependent equation of state of the form, 



P(r iif P , 6) - o 

D • J U 



such that the relationship between density, p, and temperature, 
0, yields the solution for the total hydrostatic stress T n at 
equilibrium (Figure 5.3b). This equation is called a 'kinetic' 
equation of state, and is by definition, a constitutive equation, 
At equilibrium conditions, 

T 

3 n * 5.31 

Stress is a linearly additive quantity, so the net applied 
stress, a n , may be written as the difference between the total 
hydrostatic stress, T- f , and the thermodynamic stress, p, 
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Kinetic Ecruations of State for The Gases 

Kinetic equations of state for gases are relatively well 
defined, at least at equilibrium conditions. For a perfect gas 
the thermodynamic pressure p (or T n ) can be written as, 

v 5.33 

where R is the gas law constant and V = 1/p. For an ideal gas 



m 

y equation 5.33 becomes, 



PV-pV-nRQ ^ 

5 . 34 



S3 

lu 

ltn where n is molecular weight and V is volume. If the gas is real 

P 

H then the equation of state may be defined as, 

P V - p V - Z(P,Q) n R 9 

5.35 

Z(P,0) is a correction factor from the ideal gas condition, a 
function of both pressure and temperature. The mass in a 
15 reference volume of gas does not change with variations in that 
volume, so equations 5.33 and 5.35 may be defined in terms of 
density. Letting V = 1/p we have, 
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P - p - p Z(P,Q) n R 6 • 

5.36 

and so on. For isothermal conditions, or the reversible adiabatic 
condition, the equations of state become independent of 
temperature. This is the 'barotropic' flow condition, where 
equations of state are defined solely in terms of mechanical 
variables. 

Now, the above equations of state are defined in a time 
independent fashion. This does not mean that there are no time 
dependent effects, only that these functions reflect known 
relationships between individual variables at equilibrium. 



Hydrostatic Constitutive Functions 

A hydrostatic stress applied to a volume of a simple fluid, 
regardless of the flow condition, will result in a stress limited 
change in material density. Let some fluid, a gas, be contained 
at some constant temperature 0, pressure of magnitude P c , density 
p Q , and be governed by the kinetic equation of state 5.33, 



At some other time, some other pressure P 2 , and some other 
density p 2 , let the constitutive equation be, 
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P 2 - p 2 R 6 - p 2 q 

D.Jo 



where q = R0 is constant at constant temperature* Define P 2 in 
terms of differential changes in pressure, i.e., 



P 2 - P c + dP 
2 ° 5.39 



Equation 5.24 yields the density expression and equation 5.38 may 
then be rewritten as, 



P Q + dP - (p 0 + dp) q 

Subtracting equation 5.37 yields, 

dP - dp q 

From equation 5.27 we then obtain, 



5.40 



5.41 



dP - dp q - -p p ? , , , ^ - -p 0 dVqr 



dV 

(l + dV) 5#42 



dP is the net applied hydrostatic stress, a fi , as defined in 
equation 5.32 and from which we may write the equivalence, 



- 3 dp - o n + o 22 + o 33 - 30^ 

J • 4 J 
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Now, for a perfect gas, a function K(p,0) may be defined as, 

*(..) - PoQ - Po*° . AA 

5 . 44 

The expression for small volumetric changes is given by equation 
5.23. The hydrostatic strain, € n , is equal to the mean normal 
strain, e v /3, and dV = e v * = 3€ if . From equations 5.43 and 
5.23 equation 5.42 produces, 

a n ~ - K(. .) e v - -3 .) 



which is identical in form to the hydrostatic constitutive 
function for isotropic linear elastic solids, or isotropic linear 
viscoelastic solids evaluated at t->*. The (-) indicates 
compression. For solids the coefficient function K(».) is the 
bulk modulus function, and for liquids K(..) is the bulk 
compressibility function. Equation 5.45 may also be written as, 



Now K(..) is not constant over large ranges in volumetric 
deformation. For real gases, this non-linearity is evident from 
equation 5.35. However, any non-linear function can be 
approximated by a linear function so long as the domain of 
interest is confined to an appropriately narrow range. Equations 
5.45 and 5.46 will be acceptable approximations so long as 



155 



interest is confined to small, isothermal changes of some 
reference volume. For simple fluids subjected to isothermal 
hydrostatic volumetric changes equation 5.46 also represents the 
barotropic equation of state. For small volumetric strains it is 
also a linear 'elastic' equation, where the magnitude of any 
volumetric change is a linear function of some increase or 
decrease in net applied stress. 

5.7 The Kinetic Equation of State as a Viscoelastic Function 

The Frictionless Elastic Fluid 

The simplest possible fluid is the perfectly 'elastic' 
fluid. By definition this fluid is frictionless and has no 
viscosity. The barotropic equation of state is defined in terms 
of pressure and density (volume) only, 



P ' P(P) 5.47 



The stress function for both the deviatoric and hydrostatic 
stress fields is simply, 



which states that at equilibrium the state of stress is always 
defined in terms of the thermodynamic pressure, equal to the 
hydrostatic pressure, regardless of the state of motion. No 
viscous effects are allowed, so there are no energy losses during 
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deformation, either in dilation (volumetric) or in distortion 
(shear) . While shear deformation still occurs, the only 
resistance to any deformation is in dilation. This fluid is then 
perfectly 'elastic' in the sense that volumetric deformation is 
constrained for any applied stress and completely recoverable 
under isothermal and/or adiabatic conditions. No such fluid 
exists. 

The Real Viscous Fluids 

The interaction between bonded atoms under an applied stress 
provides for the elastic impulse and the delayed elasticity 
responses in solids. This interaction must occur through 
interaction between electromagnetic fields. However, such 
interference and interaction need not require bonding. An example 
of unbonded molecular/atomic interaction is to place a heavy 
solid object on top another solid supporting object. The atoms in 
each object at the interface need not be bonded for the lower 
object to provide the counteracting support. Likewise, the atoms 
in each object at the interface are not physically touching. 
There is some spacing between the atoms that is determined and 
maintained by the electromagnetic repulsive forces acting between 
atoms. These load induced forces are then transmitted throughout 
the objects through the interaction of the electromagnetic fields 
of each bonded atom. 

The interaction between unbonded polyatomic molecules in 
fluids occurs similarly, but the unbonded molecules are usually 
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in a state of random motion and orientation. Consider the case of 
one fluid molecule colliding with another fluid molecule, 
resulting in the transfer of momentum from one moving molecule to 
the other. Momentum transfer comes about as a result of some of 
the kinetic energy being converted to elastic potential energy 
through the interaction of the electromagnetic fields between the 
atoms of the two molecules, including bond distortion* That 
stored elastic potential energy in the distorted molecule is then 
released through a repulsive (elastic) rebound reaction and, in 
combination with the electrostatic repulsive forces, etc. , acting 
between the molecules causes each molecule to assume a trajectory 
that causes collisions with other molecules which in turn rebound 
and collide with other molecules, and so on. 

Recall that the elastic impulse response is mandatory in 
solids because stress transmission cannot possibly be 
instantaneous. The same holds true for fluid molecules as well 
because the manner in which bonded atoms in the unbonded 
molecules of a fluid connect and interact is the same. In fact, 
cooling a fluid results in the bonding of the individual fluid 
molecules to form solid crystals. The chemical makeup of the 
molecules does not change with phase changes, and the atoms of 
each molecule remain bonded. Interaction between molecules of the 
solid phase or the fluid phase must therefore be described by the 
same form of constitutive function. The values of the parameters 
in this function may vary, but the basic form of the function 
will not change. So, the interaction between atoms and the stress 
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induced straining of the bonds of a fluid molecule will be 
defined in terms of a viscoelastic constitutive function with two 
basic component responses: an impulse response, and a delayed 
elasticity response. Because fluid volumetric deformation is 
stress limited, constrained by the thermodynamic and 
electromagnetic forces, the generalized solution for a 
constitutive function may not, by definition, be that of a 
'fluid'. A function analogous to that of a viscoelastic solid 
must describe volumetric deformation. 

The Generalized Hydrostatic Constitutive Function 

Equations of state are determined from pressure and 
density /volume measurements, at some temperature, all of which 
are measured at an equilibrium condition. This is analogous to 
deriving constitutive equations for viscoelastic solids by 
measuring the ultimate volumetric strain associated with some 
load and ignoring time dependent effects (e.g., Hooke's Law of 
linear elasticity) . The most generalized form of any linear 
viscoelastic hydrostatic function evaluated at t->» may be 
written as, 

ait-**).. - [f(t-<») ijkl ] e oikk) 6^ 

5 .49 

The * 0(kk) *6 kl are volumetric strains, and the a(t->«o).. are the 
hydrostatic stresses, the mean pressure, P. 

For solids we demand coupling between the initially assumed 
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independent f(t->») i j kl functions through bonding. That is not the 
case for the simple fluids. The deviatoric constitutive behavior 
of fluids is markedly different that the hydrostatic constitutive 
behavior, a consequence of the lack of bonding, or any other 
coupling forces of significance acting between molecules. Because 
there is no bonding we may not even demand coupling between the 
deviatoric constitutive functions. Unless it can be shown or 
argued that a fluid is isotropic/homogeneous, or behaves as such, 
we must admit the possibility that deformation under shear will 
be described by separate deviatoric functions, each one of which 
is associated with one of the shear stresses. An example of a 
non-isotropic fluid would be a polymeric fluid with the long 
chain molecules aligned preferentially in one direction. 

The lack of bonding does not however, preclude coupling in 
the hydrostatic f(t->*) ijkl functions. Coupling does occur and the 
coupling mechanism is pressure, the result of the thermodynamic 
and electromagnetic forces acting between the unbonded molecules. 
When coupling may be demanded equation 5.49 always reduces to, 

o(t-«) n - [d ijkl ] e oikk) b kI f(t-co) 

o • dU 

The previously discussed equations of state normally demand 
material isotropy, because the molecules of a liquid or gas are 
naturally in a state of random motion and random orientation. 
Expanding equation 5.50 for isotropy yields, 
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of isotropy. Given that a n = = a^, and £ n = = €33, the 
solution for equation 5.51 reduces to,. 



o u (t~*>) - f(t-«) <d xl + 2d 12 ) 6^ 



5,52 



or from equation 5.46, 

P(t-oo) - K(t-°o) e v 



5.53 



where the solution for K(t->«) is, 



5.54 



which is identical in form to the solution for the bulk modulus 
function of viscoelastic solids. 

Now consider a perfect gas. From the definition of equation 
5.44 equation 5.54 implies the following, 

5.55 



If the gas were 'real' then from equation 5.35 the equivalence 
would be, 
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/C( t-«, P, 0) - p 0 ZlP,6)n R 6 



5.56 



a non-linear function because Z(P,0) is non-linear. 

For isothermal conditions equations 5.55 and 5.56 represent 
barotropic equations of state and K(. w t->«) is, in general, non- 
linear. However, if volumetric deformations about some reference 
volume are small, K(..,t->«) may be approximated as a linear 
function. This implies that elastic components in the function 
K(,.,t->«), which provide the limit to volumetric deformation, 
may also be approximated by a set of linear constants. 

Oscillatory Volumetric Deformation 

Consider the frictionless and perfectly elastic fluid, where 
some reference volume is subjected to some small oscillatory 
deformation. Let the resulting change in density in time be 
defined by, 

p - p 0 cos (<*> t) 

5.57 

b> is frequency, in cycles/ sec. t is time. Define the equation of 
state, at isothermal conditions, to be that of a perfect gas, 
given by equation 5.33. The mean stress or pressure, P, is then, 

Pit, 6) - p 0 cos(a>t) f P c - p c R6 

5 • 58 
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Let volumetric strain be e v = e 0 *cos (u>t) . From equation 5.46 we 
obtain, 

P(t,6) - P 0 cos(w t) - e 0 K(..) cos(o t) 

5 • 59 

where K(..) is the linear time independent bulk compressibility 
function. 

If the linear bulk compressibility function, K(..), is not 
perfectly elastic and the fluid has energy losses associated with 
volumetric deformation, then K(..) is a viscoelastic function. 
The solution for net hydrostatic stress, a fi (t), as a function of 
the oscillatory variation in volumetric strain then becomes a 
function of the hereditary integral. 



Frequency Dependent Hereditary Integrals 

Recall the linear hereditary integrals of the solids, 



— ^ 5.60 



and 

When loading and deformation are defined in terms of frequency 
dependent sine and cosine functions the hereditary integrals 
reduce to, 
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e £ (t) - m u f* B Oj ( t) + n i;j C* (w j) o 5 ( t) 



5.62 



and 



o t (t) - a ij f s e j {t) + Cio^ejit) 



5.63 



where C*(a>) is the complex compliance function, C(u>) is the 
5 complex modulus function. Both are non-linear complex functions 
of frequency defined in terms of two component functions , 

13 

:Q C*(g>) - C* x (g>) + i C* 2 (o), C(Ci)) - qCto) + i C 2 (g>) 
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The complex stress function of equation 5.62 is defined as, 



o^(t) - o oij) [cos(o)j t) + i sinfu., t) ] - a oU) e J 



5.65 



li& and the complex strain function of equation 5.63 is defined as, 



e^(t) - e oU) [cosfwj t) + i sin(Wj t) ] - e oU) e J 



5.66 



For loading defined in terms of cosine functions, the real stress 
functions are aj.(t) = a 0(j) »cos (wjt) , and the resulting strain 
solution to equation 5.62 is, 
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e t (t) - t ] f* B a oU) cos(u>jt) + 

5.67 

[ ) o oij) {C* 1 ( Wjlcosfo j t) - C* 2 (<*) J ) sin(w jt)} 

The strain solution associated with the real loading functions 
aj(t) = a^j/sinCujt) is, 

e i (t) - [ m tj )f* B a oij) sin (<* jt) + 

5.68 

[ ] o oij) {C* 2 (<* j) sin{u> jt) + C^to^coslG)^)} 

Consider an isotropic material subjected to a static 
hydrostatic load where the hereditary integral function 5.60 
degenerates to, 



e n {C) - — ° 0 <ii) 



5.69 



and the bulk modulus function, K*(t), is defined as, 
K*(t) - 3(1 + 2m 12 ) f* B + 3(1 + 2n 12 ) f* D (t) 

5.70 

- K* B + K* D (t) 



From the hereditary integral functions the frequency dependent 
solution for periodic, time dependent hydrostatic loading and 
deformation take the form of, 



&u\t) - — °o(ii) e 3 — Oiiic; 



5.71 



where the solution for K*(w) is, 
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K*(<*) - 3(1 + 21^3) f* B + 3(1 + 2n 12 ) £>(w) 

5.72 

- K* B + K* D (u>) 



The volumetric strain is e v = e^, the net normal stress, a n (t) , 
and hydrostatic pressure, P(t). Following equation 5.46 equation 
5.71 becomes, 



e v (t) - -JC*(o>) P Q e ittt - - A>(o>) P(t) 



For the relaxation modulus we have, 

Pit) - - JC(a>) e D e ittt - - tf(u) e„( t) 



5.74 



where K(co) is, 



a-(w) 



3(1 + 2a 12 ) 



C(o>) 
3(1 + 22> 12 ) 



5.75 



K E is obviously the classic linear elastic bulk modulus function. 



5.8 Work, Energy Loss and Propagation Characteristics 

In 'perfectly' elastic solids the work of deformation would 
be stored as elastic potential energy and completely recoverable. 
In viscoelastic solids some of the energy expended in the 
deformation process is lost, i.e., internalized and transformed 
into 'heat' through friction (electromagnetic field interaction 
effects) , and is not recoverable as mechanical elastic energy. 
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Briefly, the internalization occurs through the deformation of 
bonds and vibration between bonded atoms. Ignoring conduction, 
energy or heat loss then occurs through radiation - photon 
(electromagnetic energy) emission. Unlike the solids, gas and 
5 fluid molecules are free to translate laterally, and this kinetic 
energy of motion is transferred and internalized through 
molecular collisions. For polyatomic gas molecules rotation must 
also be considered, with rotation also producing photon emission 
(ignoring any transitions between electron energy states, as with 
10 monatomic gases) . All electromagnetic emissions are, of course, 
,q quantized. 
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The Frequency Dependent Work Functions 

Consider a differential element of a solid material being 
acted upon by some force a^A, where A is the cross sectional 



15^ area of each face of this element. The dimensions of the element 

H j 

in are defined to be dx f . In a differential time interval, dt, these 
H forces cause an increase in induced strain by an incremental 

strain factor (dCj/dtJ^dt. The work, dW, performed per 

differential unit of time, dt, is then, 

o i A dx i -^f dt - dWAdx i ^ _ 

20 1 1 dt 5.76 



For any finite time interval the solution for mechanical work W 
is the integral of this equation, 
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» - / M - / °< St dt 



5.77 



(In general, this expression is not the solution for total 
work, for there is an additional kinetic energy term, call it KB, 
that comes about from the general consideration of both the 
effects of surface tractions and body forces. Letting 2 be the 
total work expression, the most general differential form of the 
total work function has the form of, 

dS _ dKE + dw 

dt dt dt 5.78 



where equation 5.77 is the integral, in time, of the differential 
mechanical work term dW/dt. The kinetic energy of a solid body 
subjected to dynamic deformation is, 



1 r de i dei i de i de i 

~2 Ja ~dt ^dt 2 ? ^t ^t 



5.79 



where p is mass density and A is area, defined here to be unity. 
The time derivatives of strain are velocity terms. The time rate 
of change in kinetic energy is then just, 

dKE _ d 2 €j Cfej 

dt " 9 dt 2 dt 5.80 



This expression plus the time derivative of the mechanical work 
function of equation 5.77 yield the complete solution for 
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equation 5.78, 



dS dKE dW d 2 e £ de i + de, 

dt ' dt dt " p d 3 t dt ' dt 5.81 



We have no interest in KE and therefore confine interest to the 
mechanical work function W. Note that for static conditions KE = 
0, and W appears to become equivalent to the elastic strain 
energy function U. As we shall see shortly, this is not exactly 
the case. ) 

From equations 5.60 and 5.62 the solution for strain is 
composed of two time dependent components, the impulse response 
and the delayed elasticity response, i.e., 

€i (t) - e^(t) + eVt) 5 82 

Substituting this expression into equation 5.77 then yields, 
f dc^t) M r de^(t) ^ t *Vt) _ 



Solving this function requires definition of the time 
dependent expressions for stress and strain. Define the loading 
function to be, 



which produces a solution for the impulse response strain given 



169 



by the first term of equation 5.67, 

e B i (t) - [ m £j ]f* B a 0(i) cos((J i t) 



5.85 



The solution for the delayed elasticity strain response is given 
by the second term of equation 5.67, 



e D i(t) - [ ] o oiJ) {C* 1 (<*> j )cos(v> j t) - C* 2 (a> j) ain(v> jt)) 

J • oo 



The derivative of equation 5.85 with respect to time is, 

d ** d i t lt) - " 1 m u ] ***** o oW) ain(«jt) 5<87 

and the derivative of equation 5.86 with respect to time is, 
d e D (t) 

-± - - [ n ±j ] o oij) { c* 1 (co J )sin(ci) i t) 

ac 5.88 

+ C* 2 (<dj) COS ((*> jt)) 

Substituting equations 5.87 and 5.88 into equation 5.83 
yields the generalized solution, 

W - - [ m A j ] f+ B w^ 0 {i)°o(j)/ t cos (Wjt) sin(<o jt) dt 

- [ jjy ] <a > a o(J) o 0(J) C* 1 (G) :; ) cos(o i t)sin(<i) i t)dt 5.89 

- [ ] G> i o o(i) o 0ti) C* 2 (a> j) j cos (w it) cos (o> i t)dt 

where only the first line represents the classic elastic strain 
energy function, being associated with the elastic impulse 
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response. 



Energy Dissipation 

The amount of energy lost per vibration cycle , or per unit 
of time, is a measure of the relative energy dissipation 
efficiency. Let's simplify equation 5.89 for hydrostatic loading 
and material isotropy. Under these constraints, the stress 
functions are equivalent and the Wj = w, and equation 5.84 
reduces to a f (t) = a 0 *cos(wt). A logical choice for the time 
period over which to evaluate equation 5.89 is one complete 
period of vibration, defined as t = 2tt/co. Equation 5.89 now 
becomes , 



The first term of equation 5.90, call it W E , is associated 
with the impulse response, 




0 




5.90 



0 




0 




5.91 



o 



which upon integration yields, 
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W* - A[m u }f* E o* [cos 2 (<ot>] *'(" 

2 J o 

- ±[ m iJ l f **°o 2 tcos 2 (2n) - cos 2 (0)] 

which is the classic potential energy function of linear elastic 
materials. 

The second term of equation 5.90, call it W° 1# is associated 
with the delayed elasticity response function, 

W\ - - [ J2~ ] o 2 0,(0)) f co cos (<ot) sin < g> t)dt ^ _ 

i ij o j 5.93 

and upon integration yields, 

- yf^ijl o G 2 C^tco) [cos 2 (cot)] 2 Y W Q4 

- ±[ ] a 0 2 O^co) [cos 2 (2tO - cos 2 (0)] - 0 

So, this function must also represent recoverable elastic 
potential energy. 

The last term of equation 5.90, call it W° 2 , is likewise 
associated with the delayed elasticity response function and is, 

2n/u 

W D 2 - - [ ] a c 2 C* 2 (co) f co cos 2 (co t) dt 

o 5.95 

which upon integration becomes, 



5.92 
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W D 2 - - [ n u ] co o c 2 C* 2 (w) [4 + -r- sin(2(ot)] 2 Y tt 
- - t n u ] w o c 2 C* 2 {w) [-* + ^ sin(4ii)] 

The term sin(4tr) = 0, so this solution reduces to, 



which is not equal to zero, and represents the mechanical energy 
lost (converted to heat) in every deformation cycle. Note that 
the solution for W° 2 is never equal to zero, except when 
t = 0, and that the magnitude of this term will increase with 
each deformation cycle. An equivalent energy loss function, 
defined in terms of the strain, would be, 



By convention, 'dissipation' is defined as the energy lost 

per unit of time, not per deformation cycle. To determine the 

loss per unit time divide equations 5.97 and 5.98 by the time 
factor 2*tr/(o and obtain, 



and 



173 



E D - [ b l3 ]^ e 0 2 C 2 (o>) 



5.100 



Both are positive definite functions where 'Ep' denotes the 
dissipated energy per unit time. Note from equations 5.97 - 5.100 
that E D is a quantized function. 

Dissipation and Speed of Propagation of a Disturbance 

The stored potential energy in a volume of deformed 
viscoelastic material does not equal the energy expended during 
deformation because some of the input energy is converted to 
heat. Consider now an impulse disturbance that travels from atom 
to atom along the length of a bar of a solid viscoelastic 
material. An initial impulse at one end of the bar imparts 
momentum, and hence kinetic energy, to the atoms at this end of 
the bar. Some of the kinetic energy involved in displacing these 
atoms and deforming the bonds with the following group of atoms 
in the bar is converted to elastic potential energy, and some is 
converted to heat as each bond is compressed or stretched. When 
the following group of atoms rebound from the deformation, the 
recoverable elastic energy transmitted through the bonds to the 
next group of atoms is less that the amount of energy expended in 
deforming the bonds. So, as the disturbance travels from atom to 
atom the magnitude of the elastic potential energy that may be 
transmitted continually decreases. This means that the amplitude 
of the disturbance must decrease as the disturbance travels along 
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the bar. 

It has been shown that the propagation velocity (speed of 
sound) of an impulse disturbance in a viscoelastic solid is 
solely a function of the impulse response. In linear materials, 
this propagation speed is independent of the magnitude of a given 
impulse or time dependent disturbance. However, energy loss rate 
is very much frequency dependent. The faster the oscillation (the 
higher the frequency) the greater the rate of energy loss. 
However, it is only the delayed elasticity term that is 
responsible for this loss of elastic energy. 

In fluids the propagation velocity may be derived from the 
barotropic equation of state. For small amplitude (small volume) 
disturbances the speed of sound is given by a constant that is a 
function of the change in pressure and density (volume) with 
respect to some reference condition, i.e., 

'* " # 

<*P 5.101 



where c is the compressional speed of sound. 

For a frictionless, perfectly elastic fluid this 
relationship implies that the speed of sound may always be 
determined directly and precisely from the equation of state. The 
solution for 9 o 9 then becomes, 



175 



From equations 5.41 and 5.42 equations 5.101 and 5.102 may be 
rewritten as, 



3 ■ " ■ • ^ ' >'*> 5.103 



and so on for ideal and real gases. As a frequency dependent 
disturbance travels through the perfect frictionless fluid and 
away from a source location there will be no energy loss, and 
hence no decrease in disturbance amplitude or intensity (aside 
from spreading effects) . 

Let's now examine wave propagation in 'real' molecular 
fluids. Let the disturbance be small enough such that the 
hydrostatic constitutive function can be approximated by a 
barotropic equation of state. No polyatomic fluid molecule will 
be 'perfectly' elastic so the equation of state must be given by 
some sort of viscoelastic function. 

In order to propagate a disturbance in a viscoelastic 
material there must be an impulse response. For small 
disturbances the speed of propagation should be relatively 
constant. However, because the constitutive response must be 
viscoelastic, there should be energy loss associated with a 
delayed elasticity response component, and the energy loss rate 
should be a non-linear function of frequency. Furthermore, as the 
disturbance propagates there should be a reduction in the 
amplitude/ intensity of the disturbance. 

Table 5.1 gives experimental data for the speed of sound in 
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dry air at several frequencies. As can be seen, the speed of 
sound is virtually constant, and is definitely frequency 
independent. These data also show that for small disturbances the 



Table 5.1 



Wave Propagation Characteristics 
of Still Dry Air at 68° F (2 0° C) 



Frequency Velocity Absorption 

(Hz) (k-ft/sec) (Db/sec) 



0.174 
0.368 
0.433 
0.573 
0.631 
0.672 
0.780 
0.970 
1.423 
3.039 
9.032 
12.306 
17.923 
25.901 
91.759 
354.700 

Note: Dry denotes the relative humidity = 0%. 
k-ft/sec denotes 1000 's ft/ sec. Hz is 
frequency in cycles/sec. Db is decibels. 

1979 CRC Handbook of Physics and Chemistry 
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5 impulse response can readily, and quite accurately, be 

approximated by a linear elastic impulse term to the governing 
constitutive equation, the barotropic equation of state. These 
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data also show that there is indeed energy loss with propagation, 
and that the energy loss is very much a function of frequency. 

So, if the state equation is described in terms of a 
viscoelastic function it is possible to account for a host of 
observed behavioral characteristics in polyatomic fluids 
including: limited volumetric deformation for any net stress, 
volumetric rebound after stress release, wave propagation 
velocities that are dependent upon an elastic impulse response, 
the decay in wave amplitudes/ intensities with increasing distance 
from source, and energy loss as a function of volumetric 
deformation and frequency. 

From these data, and the requirement for uncoupling between 
the hydrostatic and deviatoric response, we conclude that the 
complete form of the constitutive equation for the simple fluids 
must be written in terms of two independent and linearly additive 
responses , i.e., 



or for the Newtonian fluids, 



5.105 



where the pressure function 'P' is defined by a viscoelastic 
equation of state. 
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Extensions to the Atomic Level 

Consider now, for example, a chemically bonded diatomic gas 
molecule oscillating over a constant range in amplitude 
£ = c o «cos(o>t) and frequency cj. The viscoelastic work functions 
5,97 and 5.98 are energy loss terms that are functions of 
frequency, o>. This implies that if energy loss (which is 
quantized in terms of the classic vibration frequency w) is in 
the form of photon emission, then the energy of the emitted 
photon, always defined in terms of frequency, should be a 
function of this oscillation frequency, w. The energy of the 
emitted photon would then be, 



photo* 5.106 



where h is Plank's constant. This result is in accord with 
observation. However, quantum theory follows the classical theory 
of elasticity and assumes that the harmonic oscillation of bonded 
atoms is 'perfectly' elastic. This assumption precludes any 
prediction and requirement of energy loss, or the functional 
relationship of that energy loss with the frequency of molecular 
vibration. 

Interestingly, because of the assumption of 'perfect' 
elasticity for harmonic oscillation between bonded atoms, gases 
are often assumed to be 'perfectly' elastic, or very nearly so, 
although this clearly cannot be the case. (In texts physicists 
often assume 'perfect' elastic behavior while engineers generally 
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assume /perfectly' viscous behavior, i.e., the Navier-Poisson 
constitutive function for the Navier-Stokes equations.) When a 
collision occurs between two molecules of a polyatomic gas a 
'perfectly' elastic reaction can only occur if there is no net 
decrease in the sum of the kinetic energies. Assuming no change 
in rotational energies, which are quantized, a collision between 
molecules may cause a change in the non-quantized amplitudes of 
interatomic vibration. Any change in amplitude produces a change 
in vibration energy even though the vibration frequency, which is 
quantized, may remain unchanged. So, some of the translational 
kinetic energy is internalized within the molecule and eventually 
dissipated as 'heat', i.e., photon emission. A perfectly elastic 
reaction may occur only for a relatively unique matching of the 
various ienergy states of molecules at the moment of collision, a 
result that is probably very rare for even the simple diatomic 
gases. 

We may further observe that if polyatomic gases were indeed 
'perfectly' elastic then the equations of state should be 
accurate predictors of the speed of sound, for a broad range of 
temperature and pressure conditions, and energy losses should be 
negligible as a disturbance propagates. There would also be no 
reason to presume that when liquified or solidified these 
polyatomic compounds should then not also be 'perfectly' elastic, 
and the bulk compressibility or bulk modulus function an accurate 
predictor of the speed of sound. However, this is not the case. 
For heavier and more complex polyatomic gases, such as the 
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hydrocarbon gases methane, propane, etc., the equations of state 
are notoriously inaccurate predictors for the speed of sound. 

'Perfectly' elastic interactions between polyatomic 
liquid/gas molecules can only be rare occurrences. However, the 
same can not be said or interactions between monatomic gas atoms 
and/or atoms and subatomic particles (e.g., electrons). Perfectly 
elastic interactions are known to readily occur if the collision 
energy between atoms, or atoms and subatomic particles, is not 
sufficient to excite the electrons of a monatomic atom to higher 
energy states, which are always quantized. Classic experimental 
examples are the interactions between electrons and hydrogen 
atoms, and between electrons and mercury atoms in a vapor. 

5.9 Viscoelastic and Viscoplastic Fluids 

Viscoelastic Fluids 

Not all fluids are simple or Newtonian. Some of these 
materials, which have no bonding between molecules, and are 
therefore thought of as fluids, do have the ability to support 
some level of deviatoric stress. Coupling between molecules may 
be mechanical or electrostatic, or both, but these mechanisms 
provide for an interaction that supports a deviatoric stress. 
Examples of such complex fluids are the liquid uncross linked 
polymers used in the manufacture of rubbers, plastics, etc. 
Figure 5.4a illustrates the observed constitutive response under 
a unit stress. The molecules of such polymers are long chain 
hydrocarbon molecules. These fluids are often conceptualized as 
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being analogous to a squirming, entangled mass of worms, where 
the physical entanglement of the molecules provides mechanical 
coupling without bonding. 

Now, a pure shear stress may always be converted into a 
normal deviatoric stress, and we see from Figure 5.4a that the 
elastic impulse strain and the delayed elasticity strain for such 
a stress is recoverable. So, a constitutive function for such a 
material would essentially be defined like the constitutive 
response of a solid, except for the unconstrained creep response. 
The elastic impulse response and the delayed elasticity response 
are defined and coupled for both hydrostatic loading and shear 
loading, but the unconstrained creep response can only be defined 
in shear. 

Viscoplastic Fluids 

Another class of materials, the so-called viscoplastic 
fluids, again have no bonding between molecules, but mechanical 
and/or electrostatic forces again act to couple the unbonded 
molecules and prevent unconstrained deformation until some 
minimum level of shear stress has been obtained. Examples of such 
materials are drilling muds. These highly useful materials are 
often simple concoctions of water, starch (a long chain polymer) , 
and bentonitic clays, and are used extensively in the drilling of 
oil, gas and water wells. 

The simplest of these materials is the idealized Bingham 
plastic. Figure 5.4b illustrates the relationship between applied 
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stress and a linear rate of deformation for this material. Of 
course, in general, the relationship between stress and rate of 
deformation need not be strictly linear. These materials behave 
as solids until a critical value of shear stress is reached, then 
5 behave as a fluids. The idealization of Figure 5.4b may 

illustrate the deformation properties under shear, but it does 
not account for material behavior before flow is initiated. Prior 
to flow the response must be that of a solid, specifically a 
viscoelastic solid. As such, the functions describing the 'solid' 
10 phase of behavior must have two component responses, an impulse 
J response, and a delayed elasticity response. 

:~ Other types of non-solid materials that show this same type 

^ of behavior are granulated solids, such as flour, sugar, and 

*S common beach sand. These materials may be piled into free 

ID 

1*5 standing mounds that clearly support some minimum level of shear 

!□ 

'""4 stress. The coupling mechanism between particles is 'frictional' 

ly 

|n resistance - surface roughness coupled with granular angularity 

Q 

U and lithostatic (overburden) pressure (marbles, ball bearings, 

etc., will not form free standing mounds). Until the 'frictional' 

20 resistance between particles is overcome the material behaves as 
a solid, capable of supporting and transmitting a shear stress. 
Once this 'frictional' resistance has been exceeded the particles 
of the material begin to flow in a fashion analogous to fluids. A 
mass of smooth spherical particles, like marbles, move and behave 

25 in a fashion analogous to the simple fluids. 
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5.10 Inelastic (Plastic) Deformation of Solids 

Plastic or inelastic deformation covers a host of behaviors, 
some of which are physically destructive and occur when the bonds 
in a solid are broken, crystals are cleaved, and separation 
planes form internally between continuous masses of material. Not 
all inelastic deformation in solids is destructive and not all is 
non-recoverable. Some materials, which show a permanent set as a 
result of very large strain, may return to their to the pre- 
deformation shape when heated. 

So long as the deformation process in solids is not 
physically destructive, and molecules remain bonded, we may 
demand coupling in the associated strain fields. Without 
specifying the controls over any recoverability in the inelastic 
component, a generalized linear constitutive function for such a 
material would take the form of, 

e;(t) - {irtij] £** + [n tj \ f* D (t) + [ o u ] f* p (t, ..) 

Likewise for the non-linear equations. Each function represents a 
physically distinct and measurable response. 

Inelastic deformation may result in more than the simple 
crystal lattice and bond deformation. It may include the breaking 
and reforming of bonds, producing extremely large strains. An 
example would be high temperature molding of metals. However, 
large deformations are analogous to fluid flow, and can only be 
accommodated deviatorically (e.g., the Bingham model). 
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SECTION 6 

APPLICATION EXAMPLES - ELEMENTARY BEAM THEORY & MATERIAL 
APPROXIMATIONS 

To illustrate how constitutive functions are employed, let's 
review the development of elementary beam theory. Originally 
developed using the equations of Hookean elasticity, it may now 
be extended to include the viscoelastic materials as well. 

6.1 Shear Forces and Bending Moments 

Consider any engineered structure, a house, bridge, ship, 
etc. Many of the load carrying members of these structures are 
'slender members', i.e., a structural member where the length 
greatly exceeds the lateral dimensions. A 'slender member' may be 
a beam, a shaft, a rod, a strut, and so on. These structural 
members may be bent, pulled, and twisted. In engineering design 
we are interested in the internal forces, moments, stresses, and 
the deflections and strains of the members that result from 
loading and deformation. 

When a slender member is subjected to transverse loading, we 
say that it acts as a 'beam'. Examples of slender members that 
act as 'beams' are joists and girders, leaf springs, even the 
wings of airplanes act as beams (albeit complicated beams) . 

The following discussion pertains to viscoelastic 'beams'. 
The purpose is to illustrate the method of application of the 
linear viscoelastic constitutive equations in engineering design. 
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For the sake of brevity, interest is restricted to beams with 
cross-sectional areas that are symmetric with respect to a cross- 
sectional plane located along the length and through the middle 
of the beam. Material properties are constrained to be linear, 
homogenous and isotropic, and only static loading conditions will 
be considered. 

Differential Equations 

Three forces and three moments act about a cross-sectional 
surface within a transversely loaded rectangular beam (Figure 
6,1). F w is the axial force which elongates the beam. and P u 
are the forces that shear one part of the beam relative to 
another. is the twisting moment acting about the x axis. 
and are the moments that bend the beam in the z and y 
directions. 

Figure 6.2 illustrates a differential element located 
somewhere along the beam (Figure 6.2a), and those forces and 
moments acting about that element (Figure 6.2b and 6.2c). The 
distributed transverse loading function is p(x). The magnitude of 
the total force acting along this element is p*Ax which acts 
through some point 0, where p is the average load intensity. 

Static equilibrium demands that the sum of the forces and 
the moments equal zero, yielding the following solution for the 
sum of the forces, ZF , acting in the y direction, 
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£ Fy " 0 - (Fxy + AF^) " Fjcy * P Ax 

D.J. 



and the solution for the sum of the moments, 2M 0 , acting about 
the z axis and in the x direction, located through point 0, 

£X - 0 - (M„ + A#f„) ♦ (F^ + AF^) ♦ AF^ - 



5 These equations reduce to, 

Ax 6.3 

=..0 
in 

and 



s 

M 

Sri 



+ F„ - — - 0 



Ax ^ 2 6.4 



As Ax->0 these functions reduce to, 



^ * P - o 



10 cbr 6.5 



and 



^ ♦ F„ - 0 



<i* ^ 6.6 



where the shear force AF^ -> 0 as Ax -> 0. These two equations 
are the governing differential equations relating load p(x) to 
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the shear force F^x) and the bending moment M^tx) . Integrating 
from x t to x 2 produces the solutions for the shear force and the 
bending moment , 

*i 6.7 

and 

W^Uj) - Af„<jq) + | F^dx - 0 

x a 6.8 



6.2 Viscoelastic Constitutive Equations 

We are interested in beams composed of all possible 
isotropic linear viscoelastic materials. Loading conditions are 
static, so the generalized hereditary integral functions 
degenerate into the functions for static loading, 

e^t) - [w^ f* E + n u f* D (tj) ] a j 

6.9 

and 

3 3 3 3 6.10 

which may also be written as, 
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- [ k^it) ] {f* B + £* D (t i )\o j 

- I k^it) ]f*{t 1 )o J 



6.11 



and 



where 



o.i(t) - t luU) ] If* * fltj] ej 

6. 12 

- [ l^U) ]f(t l ) ej 



*« lt) /.'♦/.*<e 1 > 6.13 



and likewise for l^Ct). Because only isotropic materials are 
considered, the matrices m-, n iJf etc., reduce to, 

1 m 12 m 12 0 0 

1 0 0 

.... 1 0 0 

H 0 

H 

and so on. Recall that /i = (1 - m 12 )/2. 



o 

0 
0 
0 
0 

n 



6.14 



6.3 Bending Induced Stresses and Strains 

Shearing forces and bending moments act on any cross-section 
surface (Figure 6.1). Let's now determine the distribution of 
internal stresses induced by pure bending, the bending moments 
M u (x) and M (x) . Deformation and strain are time dependent, 
reaching a maximum value at appropriately large time. So, 
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interested is initially restricted to solutions for f*(t) and 
f(t) at very large time, call it t->«, where transient 
deformation has ceased. 

Curvature Concepts 

When loaded a straight beam will deform into a curved shape. 
'Curvature' is defined as the rate of change in the slope of the 
beam curvature. Curvature may be thought of as the shape of a 
line or fiber imbedded within the originally undeformed straight 
beam and oriented parallel to the longitudinal x axis. Figures 
6.3a and 6.3c illustrate the curvature of the line segment XY of 
length AS. The change in the slope of the angle 0 between point X 
and point Y is A<p. When the change in slope angle is small, the 
small angle approximation yields the solution for arc length AS 
as R # A0. Curvature of a line segment, d$/dS, is then, 

d<$> _ lim A<t> _ lim _1 m JL 

ds " AS-o AS " AS-o r ~ r 6.15 

where R is the 'radius of curvature' at point X. 

Consider the geometry of deformation of a beam illustrated 
in Figures 6.3a and 6.3b. The lines A, B, and C represent three 
equally spaced cross-sections through the undeformed beam. A', 
B', and C represent the same cross-sections in the deformed 
beam. During bending the two sections of the beam, defined by the 
3 cross-sections (A, B, and C) , are deformed. Because of symmetry 
in the bending moments about the x-y plane of symmetry, 
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deformation should be symmetric about that plane. Let the cross- 
sections be placed such that the two beam elements defined by the 
three cross-sections are identical. By demanding that the 
transmitted moment throughout the beam be constant, it is 
reasonable to expect that the deformed shapes will likewise be 
identical. However, changes in shape of the two elements will be 
accompanied by a corresponding changes in shape of the three 
cross-sections. Any such changes would require shape changes in 
the surfaces defining the sides of contact between the two 
adjacent elements. Furthermore, moment symmetry requires symmetry 
in deformation, and therefore mirror image shape changes in the 
contact surfaces, a physical impossibility. Hence, all cross- 
sections must remain plane and normal to the plane of symmetry, 
and each element deforms identically, due to the constant M X2 
along the length of the beam. Extensions into space of the 
initially parallel cross-sections will then all have a common 
intersection point, and the beam bends into the arc of a circle 
centered about this point. 

The above arguments have not ruled out deformation of a 
cross-sectional plane within it own plane, which does occur. The 
only restriction imposed on such deformation is that it must be 
symmetric with respect to the plane of symmetry of the beam. 

Deformation Geometry and Strain 

Figure 6.4a illustrates a segment from a beam before 
deformation. While any cross-section through the beam remains 
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plane during deformation, any straight line or fiber parallel to 
the x axis is bent into an arc. Some will elongate while others 
shorten. There is however, some location within the beam, as yet 
unknown, where these lines remain unchanged in length. This line 
is called the 'neutral axis', and our reference coordinate system 
is set up so the x axis is coincident with this 'neutral' axial 
line. By convention, the x-y plane is the 'plane of symmetry', 
and the x-z plane is the 'neutral surface'. 

We are developing a small displacement theory. So, while 
there may be deformation of each cross-section within the plane 
of the cross-section, deformation is constrained to be 
sufficiently small so the location of a point in the undeformed 
cross-section is an adequate approximation for the location of 
that point in the deformed cross-section. From Figure 6.4 observe 
that the lines AB and CD, separated by some distance y before 
deformation, are deformed into the arcs AB' and CD' after 
deformation. Given the above constraint these lines will be 
separated by approximately the same distance, y, after 
deformation. The radius of curvature for the line AB' is R, but 
R - y for the line CD', y is measured from the neutral axis. 

Since the neutral axis is a line of zero deformation, the 
lines AB = CD = CD', and yield the strain solution, 

AB f - AB AB* - CD 7 

AB CD* 6.16 

Using the small angle approximation, sin0 * 0, and letting each 
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line represent a differential quantity, the arcs formed by the 
lines may be expressed in terms of dtp, 



CD' - R d$> , AB' - (R - y) d$ - . _ 

O.J./ 



Let dS = CD', From equation 6.16 the solution for longitudinal 
strain along the plane of symmetry becomes, 

R dS Y 6.18 

The minus sign indicates shortening above and extension below the 
neutral axis. From symmetry arguments requiring plane sections to 
remain plane we conclude that the shear strains = € JiZ = 0. 
Because the transverse shear strains = = 0, it follows 
from equations 6.9 and 6.14 that the corresponding shear stresses 

and must also equal zero. However, no definitive 
statements may yet be made about the strains e^, € zt9 and 
beyond the requirement of symmetry about the plane of symmetry. 

Additional Equilibrium Equations and Stress Solutions 

Equations 6.9 and 6.14 produce the constitutive equation for 
isotropic materials that relates the maximum longitudinal strain, 
€ x (t->») , to the appropriate stresses, i.e, 

6.19 

+ f* D (t~~) [n^ a x + n 12 (o y + a s )) 
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or after rearrangement and following equation 6.11, 



y 6.20 



where from equation 6.13, 

f* s m 12 + f* D (t^) n 12 



k 12 ( t-«) - 



f** + f**(t-«) 6.21 



and m n = n n = k t1 = 1, by definition. When the delayed elasticity 
component of the constitutive function is negligible, f*(t) 
approaches the linear elastic solution, i.e., f*(t) -> 1/E and 
k 12 -> m 12 -> -w. E is Young's modulus, and v is Poisson's ratio. 
So, equations 6.11, 6.20 and 6.21 are not only identical in form 
to the linear elastic solutions, they yield the linear elastic 
solution as f* D (t) -> 0. 

Equation 6.18 provides a method of solving for the 
longitudinal strain, € x (t->») , if y and R can be measured. 
Unfortunately, equations 6.19 and 6.20 do not provide unique 
solutions for the three normal stresses, a y , a z , and o x . 

Because we are considering pure bending only the shear 
forces P^ and equal zero. If dA is some differential area on 
the face of a cross-section through a beam, then the force acting 
on dA is simply dP x = a x »dA (Figures 6.5a and 6.5b), where 
SP X s 0 at equilibrium. This constraint yields, 
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u 



s 



A 6.22 



Symmetry in the loading force (P x ) about the plane of symmetry 
yields the equilibrium equation for the bending moments My, 



E"y - / z o x dA - o 



6.23 



Conversely , loading symmetry about the neutral surface (the x-z 
plane) has not been required. So, the corresponding equilibrium 



U equation for the bending moments M w becomes, 



ir| * ** * ** 6.24 



None of the above expressions yield any information about the 



in 

1CH transverse normal stresses <r , <r, or the transverse shear stress 

ill 

\*k However, the external surfaces of the beam are free of such 

a cs. 

It forces, due to the imposed boundary conditions. Given that the 
beam is slender, it is reasonable to assume that any such 
stresses, if they exist, are much smaller than the induced 
15 loading stress, a x , and negligible in comparison. We therefore 
impose, 



y 6.25 



throughout the interior of the beam. These conditions leave only 
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one normal stress, a x , and equations 6.19 and 6.20 reduce to, 
e,(t-») - [f** + f* D (t~<°)] o x - f*(t-«) o x 



6.26 



Combining this expression with equation 6.18 yields the solution 
for the longitudinal normal stress a x , 



JL dfc 



R f*{t^<») f*(t-~) dS 6.27 



which is identical in form to the linear elastic solution, and 
;2 will degenerate to that solution when f* D (t) is negligible. 

tfl Moments of Inertia and the Flexural Rigidity Term 

y 

IB Substitution of equation 6.27 into equation 6.22 produces, 

if? 

O TF X - f o x dA - f p . y . dA - _ . ). , f ydA - 0 

IOlJ Ja fif»(t-») R f * ( t-«) J a 6.28 

IU 

in 

D 

1,4 The final integral term is the first moment of the cross-section 

p. ■ 

area. R and f*(t) do not equal zero for any real material or any 
bent beam, so the integral term must equal zero, a result that 
demands the neutral axis pass through the centroid of the cross- 
15 section of the beam. 

Substituting equation 6.27 into equation 6.24 yields, 



J>. - [ A yo x dA . R f a y><* - M bt 



6.29 
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where (M^) represents the solution for the applied bending 
moment about the z axis. The integral function on the right hand 
side is the second moment of the beam cross-sectional area, the 
'moment of inertia' 1^, defined by, 



6.30 



Redesignating f*(t->«)/i yy as the 'bending modulus' or 'flexural 
rigidity' of the beam, D(t^><»), yields, 

I yy 6.31 

When f* D (t) -> 0, then f*(t) -> 1/E, the linear elastic solution, 
and then D(t) -> 1/EIyy, the standard definition of D for linear 
elasticity. 

Final Stress and Strain Solutions 

Substituting equation 6.31 into equation 6.29 yields the 
solution for 1/R, 



1 

* ~ DZ ~- ' 6.32 



4 - M bz Z)(t-«) 



Equations 6.18 and equation 6.32 then produce the final solution 
for the longitudinal strain, € x (t->») , in terms of M^, 
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e^t-o) - -M bt D(t—>)y 

6.33 



The solution for longitudinal stress, a x , is obtained from 
equation 6,27, 

M bz y ejt-o) 



Given that a y = a z = 0, equation 6.21 and equations 6.9 and 6.11 
produce the following solution for the transverse normal strains, 
€ y (t->») and e z (t->oo), 

y 12 x 6 35 

Substituting equation 6.34 into this solution yields the final 
solution for the transverse strains defined in terms of the 
applied bending moment, 

€ x (t-») - e y (t-«) - - k 12 (t-oo) M bz D(fc-«) y- jc 12 (t-«)e x (t-«) 

6.36 

The strain = 0 due to the stress conditions imposed by 
equation 6.25. 



6.4 Solutions for Deflections 

Bending Moments and The Moment-Curvature Function 

When a symmetric, linear viscoelastic beam of isotropic 
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material properties is subjected to pure bending the curvature of 
the neutral axis is defined by, 

* - "I - "-^^ te 6.37 



If a beam is conceptualized as a connected series of differential 
elements then equation 6.37 describes the curvature of the 
neutral axis in each of those elements. Extending this function 
to a generalized loading condition we make the approximation that 
the internal shear forces, which necessarily accompany any 
variation in bending moments, do not make any significant 
contribution to bending. Given this approximation any variations 
in beam curvature may be determined from equation 6.37. 

To determine the curvature of the bent neutral axis it is 
necessary to deduce its displacement, defined as w(x) , from a 
static x axis (Figure 6.6). The slope of the neutral axis may be 
defined as, 

- tan* 

6.38 

Differentiating with respect to the arc length S yields, 

d , dw(x) * d 
-dS { -dbT ) ' dS (tan *> 6 .39 

Making use of the differential relationship, 
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d(. .) d(. .) dx 



10 

p 
Si 

IU 

in 



dS dx dS 6.40 



produces the differential function, 



d^wUldx m sec 2^d$_ 
dx 2 dS * dS 6.41 



Solving for the curvature d0/dS then yields, 



d<j> . ^MxL_^ cos 2 <t) 

dS dx 2 c?5 6.42 



o 

*Q From Figure 6.6 the solution for the cos0 function is, 



. dx 
cos* - - 



dS [1 + (dv(x)/dx) 2 ] J/3 6.43 



Substituting this expression into equation 6.42 produces, 



d<t> _ dMx) /dx 2 

d5 " [i + (dw(x)/dx) 2 ] J/3 6.44 



10 and from equation 6.37 we obtain, 

1 dV(x) /dx 2 



" Z>(t-~) [i + (dt/(x)/dx) 2 ] 2/3 6.45 



Restricting interest to small angle deflections, say on the 
order of 5 degrees or less, then tan0 = dw(x)/dx has a magnitude 
of 0.1 or less, and [dw(x)/dx] 2 has a magnitude of 0.01 or less. 



200 



13 

IU 
111 

a 

I" 



So, neglecting this squared term in equation 6,44 yields, 



d 2 w(x) m _d£ 
dx 2 dS 6.46 



which has an error of one percent or less when compared to the 
exact solution. Equation 6.45 then reduces to, 

M - 1 d 2 w(x) 
bz D(t-oo) dx 2 6.47 



and is known as the moment-curvature relation. 



*D The Load-Deflection Function 



Refer back to equations 6.5 and 6.6. Differentiate equation 
6.6 again with respect to x and substitute those results into 



10 equation 6.5. This yields, 



^ . P 

dx 3 6.48 

Combining this result with equation 6.47 produces the expression 
relating loading to displacement, 

d a ( 1 d 2 w{x) . _ 

dx 2 D(t-«) dx 2 6.49 



15 If D(t) is constant along the entire length of the beam, i.e., 
I is constant and the constitutive properties are always 
described by (t) ]*f * (t) , then equation 6.49 reduces to, 
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Equations 6.5 and 6.47 now yield the following expression for the 
shear force F in terms of displacements, 

F - 1 d*w(x) 

D(t-*>) dx 3 6.51 



6.5 Time Dependent Solutions 

The displacement function w(x) has only been expressed in 
terms of the independent variable x. However, because f*(t) is 
evaluated at t->«, it is obvious that the displacement function 
is also a function of time and that w(x) = w(x,t->«>). 

Loading conditions are constrained to be static. So, how 
does a change in the time variable from t->* to t = 0+ effect the 
previous derivations and solutions? In this case f*(t) reduces to 
f*(0+) = f* E , and k 12 (0+) = m 12 . The bending modulus function D(t) 
reduces to D(0+) = f^I^. These functions represent the impulse 
response, or similarly, the linear elastic constitutive functions 
if the delayed elastic response is negligible. The only 
difference between the two solution sets are the values for 
f*(0+) and f*(t->«), etc. 

When f*(t) is evaluated at any time t it is constant valued. 
This function is the sole control over the magnitude of the 
bending modulus D(t), likewise constant valued at time t. Recall 
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that the loading forces are time independent. Therefore any time 
dependence in the solutions for the longitudinal stress, the 
three normal strains, and the resulting displacements will be 
controlled entirely by the time dependence in f*(t) and D(t) . 

Keeping in mind these observations, and that all the applied 
forces are time independent, then from equation 6.33 the time 
dependent solution for the longitudinal strain, € x (t) , becomes, 



- -M bz D(t)y 



From equation 6.3 6 the time dependent solutions for the lateral 
transverse strains, € z (t) and € y (t) become, 



e z (t) - e y (t) - - k 12 (t) M bz D(t) y - k 12 (t)e x (t) 

6 • 53 



Likewise equations 6.34 and 6.52 yield the solution for 
longitudinal stress, a x , 

I yy 6.54 

From equation 6.50 the load-deflection equation is now, 



dx* 6.55 



Integrating with respect to x four times, and evaluating at the 
appropriate boundary conditions, will produce the time dependent 
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solution for the displacements as a function of both time t and 
location x. 

Likewise the time dependent solution for the moment- 
deflection equation given by equation 6.47 becomes, 

d*w(x.t) _ D(t) M ^ 



5 etc 2 6.56 

Integrating with respect to x twice, and evaluating at the 
appropriate boundary conditions, will yield the solution for the 
displacements as a function of time t and location x. 
O There is a caveat to these time dependent solutions which we 

10S shall now discuss and illustrate by example, 

5 : 

£ 

in 

jig 6.6 Application to an Example Problem 

IB 

Consider the problem of a cantilever beam loaded and 
supported as illustrated in Figure 6.7a, Define the beam to be 

ry 

l~ rectangular, of thickness H, width B, length L, and homogenous 

u * 

li=? and isotropic in constitutive properties defined by the linear 
function f*(t) = f* E + f* D (t). 

The Loading Function 

The loading function is p(x), where the maximum load, 
defined as p G , occurs at the supported end of the beam. The 
20 loading function, as a function of location x along the beam is, 
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Pix) - p 0 | 



6.57 



where x is measured from the unsupported end of the beam. 



The Moment of Inertia and the Bending Modulus 

The solution for the moment of inertia 1^, given by 



equation 6.30, is, 



H/2 

B H 3 



Iyy - f y 2 dA - f y 2 B dy 



m 12 6 - 58 



The solution for the bending modulus function D(t) is, 



D (t) - f * U) - 12 f *(t> 

T yy 6.59 



The Limiting Boundary Conditions 

Because the beam is imbedded in the wall at x = L, the 
displacement of the beam at x = L is zero for all time. 
Furthermore, the imbedded end of the beam is restrained from 
deflecting, so the slope of the beam at x = L is likewise zero 
for all time. The two defining boundary conditions are therefore, 

w{Ltt ) - o , dw{ 5 ' t] - o 

<** 6.60 
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Equilibrium Equations and Solutions for Shear Forces and Bending 
Moments 

Refer to Figures 6.7b and 6.7c. Summing the forces F y in the 
y direction and demanding equilibrium yields , 

Z r y - (p<* - «. - o 6 61 

where the integral represents the intensity of the load 
distributed along the length of the beam. R 0 represents the 
opposing resultant force required to maintain equilibrium. 
Substituting equation 6.57 into 6.61 and integrating the result 
produces the solution for the resultant force R 0 , 

Figure 6.7d illustrates the shear force and moment 
distribution at some location within the beam. Again, summing the 
forces F y in the y direction and demanding equilibrium yields, 



£ F y - - / P dn + F 3 



6.63 



where F w is the internal shear force. Inserting the loading 
distribution function p = p G # n/L and integrating yields the shear 
force solution as a function of location x, 
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At x = L the sum of the resultant maximum moments, M L , must 
equal zero for equilibrium to exist. The force exerted by a load 
p«dx across the moment arm L - x (Figure 6.7c) produces the 
moment p«dx»(L - x) . The sum of all these differential moments, 
and the opposing moment 1^, must equal zero. This leads to the 
integral expression, 



Y,M L - 0 - f p(L - x)dx + M bx 



Q " J, 6.65 

*o 
!=* 

Iff Substituting equation 6.57 into equation 6.65 and integrating 

W 

lffl produces the final solution for the maximum bending moment, 

i 

P -(P.|a^* " ^ 6.66 

10 

in 

Refer once again to Figure 6.7d. Following the same procedures as 
before produces the solution for the bending moment along the 
length of the beam, 

15 - - Af 1 *-"'"" . 6 . 61 

which is independent of w(x) . The stress function of equation 
6.54 is therefore also independent of w(x) , and must be in order 
to obtain a time dependent solution for w(x) . 
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The Solution for Displacements, Strains, and Stress 

Substituting equation 6,67 into equation 6.56 produces the 
differential equation, 

d2 "lx> fc > - - D {t) Po -Si 

dx* Po 6L 6#68 



5 Integrating this expression once yields, 



dx ^° 24L 1 6.69 



13 Recall the constraining boundary condition dw(x,t)/dx = 0, at 
kQ x = L. This boundary condition yields the solution for the 



hi 

m 



U 

ill 
O 



constant of integration C,, 



sB C, - D{t) p. 

1Q- 1 5 24 6.70 



Substituting equation 6.70 into equation 6.69 and integrating 
again yields, 

*<x.t) - -D(t)p ol g z *D(t)p.4 i £*q l 6>?1 



Evaluating at the boundary condition v(x,t) = 0 at x = L yields 
15 the solution for C 2 , 



C 2 - - D{t) Po ^ 



30 6.72 
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Substituting this result into equation 6.71 produces the final 
solution for displacement, 



P 



m 

to 

in 



u 

HI 



in 

o 



Substitute equation 6.59 into equation 6.73 and evaluate at 
time t->». The result is the solution for the maximum deflection 
along the length of the beam, i.e., 

W{X,t } BH> 1 10L 2 2.5 J 6 .74 



;Q Equations 6.52 and 6.67 now yield the solution for ultimate 

3 - 

1/1 longitudinal strain, 



2 P Q X 3 f* ( t-co) y 



10 x LB// 3 6.75 



For a rectangular beam the neutral axis is coincident with the 
centroid of the cross-section area, which is located at the in 
the middle of the rectangular cross-section, at H/2 and B/2. y is 
measured from the neutral axis. Hence, the maximum strains occur 
15 at y = H/2 and y = - H/2, the top and bottom of the beam 
respectively. 

The solutions for the ultimate lateral transverse strains 
are obtained from 6.53 and 6.75, 
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L B H z 



6.76 



L B H* 



which also maximize along the top and bottom of the beam. 

Last but not least, the solution for stress along 
the length of the beam is given by equations 6.54, 6.58 and 6.67, 
or from equation 6.75, 



which likewise maximizes along the top and bottom of the beam. 
This stress solution is independent of the constitutive 
properties of the material, the displacement of the beam, and is 
dependent only upon the magnitude of the bending moment. 

If stress were a function of w(x) then we could only solve 
for the impulse displacement, w(x,t=0+) and the ultimate 
displacement w(x,t=«>). A time dependent stress function would 
require the use of hereditary integrals for solution. Higher 
order theories, such as plate and shell theories, produce 
solutions where stress is very much a function of the 
displacement. 



M bm {x)y 



2 y p 0 x 3 

B H z L 



6.77 



6.7 Fiber Composite Approximations 

Reinforced fiber composite materials may be constructed by 



imbedding long fibers of one type of material in a matrix of 
another type of material. Typical examples of such types of 
materials are fiberglass and carbon fiber reinforced plastics, 
steel belted rubber tires, steel reinforced stressed and 
unstressed concrete, etc. 

It is not uncommon for the constitutive properties of the 
reinforcing fiber material to be markedly different than that of 
the surrounding matrix material. Now, both types of materials are 
always fundamentally viscoelastic. However, the disparity in the 
delayed elasticity response of one material in comparison to the 
other may be such that the time dependent response for one may 
dominate to such an extent that the time dependent response of 
the other may be ignored in defining a constitutive function for 
the composite material. Steel belted rubber tires and reinforced 
concrete would be examples of such types of composite materials. 
The delayed elasticity responses for rubber and concrete are far 
more pronounced than they are for steel. In fact, the delayed 
elasticity response of steel is so insignificant that it is 
ignored for most engineering applications. 

The purpose of imbedding fibers in a matrix of another 
material is usually to limit deformation and provide a degree of 
tensile strength that the matrix alone does not possess. The 
coupling mechanism between the imbedded fibers and the matrix is 
most often mechanical, but may be chemical as well. A non- 
reinforced mass of matrix material is, like rubbers and plastics, 
typically isotropic. If the imbedded fibers are direct ionally 
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oriented, as in a 90 degree blanket weave, they will impart a 
directionality to the constitutive response of the composite 
material (Figure 6.8). A 90 degree weave will produce orthotropy. 
Assuming linearity in the constitutive responses for both 
5 the fibers and matrix materials the governing constitutive 
response for the fibers may be written as, 

e(t) i(f) - K,u>]f**m + [r2 ijif) ]f* D (t) if} 



and for the matrix, 



1QQ If the creep compliance response of the fibers is 

10 

E negligible, then regardless of the relative volume of the fibers 



""4 
IU 

m 



with respect to the matrix material we may be able to ignore that 
response component, and the resulting constitutive function for 
the composite material would be of the form, 



5 €(t) i(c) - [m ij{c) ]f* B ic) + [n ijic) ]f* D (t) {a) 6 

If the matrix material is isotropic then a 60 degree weave 
would produce isotropy for the composite material, while a 90 
degree weave produces orthotropy. Needless to say, a similar 
scheme may be applied to non-linear materials as well. 
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6.8 Discussion 

Elementary beam theory may be extended to include linear 
viscoelastic materials now that the fundamental form of the 
viscoelastic constitutive equations are known. In this section 
5 only static loading conditions were discussed, but clearly the 

above theory may be extended to include dynamic loading. However, 
an extension for the viscoelastic materials must incorporate the 
hereditary integral functions. 

From the above discussions it is obvious that other well 
10 developed, and higher order structural theories, such as those 
O developed for plates and shells, may be extended to include the 
^0 viscoelastic materials as well. Interest need not be confined to 

Id: 

ffl structural design, but may be extended to included vibrations, 

w 

13 wave propagation, etc. In fact, the field of interest includes 

IB 

15 virtually any engineering problem requiring incorporation of 

M ... 

*g material constitutive properties into solution algorithms. 

iti 

I'fi Furthermore, now that the fundamental form of the 

I s : viscoelastic constitutive functions are known the numerical 

methods, such as finite element or finite difference mathematics, 

20 may be used to develop numerical algorithms that are incorporated 
into computer software. These algorithms in turn may then be used 
to solve scientific or engineering design problems which do not 
have closed form solutions and involve objects of very complex 
geometries, material properties, and loading. 

25 Last, but not least, we assumed that the constitutive 

properties were known. However, proper characterization and 
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quantification of the constitutive properties of a given material 
is not possible without a some knowledge of the basic form the 
constitutive function may take. This basic knowledge fundamental 
to the development of the appropriate laboratory testing methods 
necessary for accurate quantification of any given material's 
constitutive functions, the degree of material symmetry or 
isotropy, the number and values of the non-zero material matrix 
coefficients, the form of thermal expansion functions and the 
values of the associated matrix coefficients, and so on. 
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SECTION 7 

APPLICATION EXAMPLE - FLOW EQUATIONS FOR COMPRESSIBLE NEWTONIAN 
FLUIDS 

It is from the generalized equations of motion for fluids 
that the Navier-Stokes equations for compressible Newtonian fluid 
flow were developed (c. 1845) . Originally developed from the 
Navier-Poisson constitutive equations for Newtonian fluids (c. 
1831) , the Navier-Stokes equations may now be corrected to 
incorporate the proper constitutive equations. 

7.1 The Continuity Equation 

Recall that the solution for the continuity condition of 
compressible and incompressible fluid flow was based on the mass 
balance function, 

it f d(pv x ) 3(pVy) + d{pv E ) 

dt ~ dx dy dz 7.1 

which defines the time rate of change in mass density (Figure 
5.2). Differentiating the right side of equation 7.1 produces, 




7.2 



or 
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- - P (V • v) 

7.3 



where Dp/Dt is the 'substantiative' or material derivative 
defined as, 

Dt dt x dx y dy * dz 7# 4 



If a fluid is nearly incompressible, or when flow is steady state 
and density changes during flow are negligible and it follows 
that Dp/Dt » 0, This condition produces the continuity equation 
for incompressible fluid flow, 



V • v - 0 

7.5 



7.2 Force, Momentum Balance and The Equation of Motion 

Equations 7.3 and 7.5 do not completely describe fluid 
behavior because they fail to consider the forces acting on, and 
within, a given volume of fluid. Because a fluid of interest may 
be in motion, it is convenient to express this force balance 
equation in terms of the rate of accumulation of momentum. For 
any given differential volume through which fluid material may be 
moving the force and momentum balance expressions may be written 
as, 
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Rate of Rate of Rate of Sum of 

Momentum = Momentum - Momentum + Forces 

Accumulation In Out 



Momentum enters a differential volume element by two 
mechanisms. The first is convection, which is simply bulk fluid 
flow. The rate of momentum accumulation is the mass flow rate, 
5 p # v, entering the differential element multiplied by the 

velocity, v, of the moving volume of material, p is mass per unit 
volume. The second mechanism of momentum transfer is the internal 
viscous shear (deviatoric) stress effects resulting from 
f 5 molecular interaction. These viscous stress effects produce the 
10 jD velocity gradients within a volume of fluid. There are two other 

fH stress fields that also act on this differential element of 

fjj 

IS fluid. One is hydrostatic pressure, which acts equally in all 

f-fl 

e axial directions. The other is a body force induced by gravity 

S3 

Ci and which acts about the mass centroid of the differential volume 

W 

15 ifi (Figure 5.2) . 

The Summation of Forces 

Let's now sum the forces acting in the x direction of flow 
through the differential fluid element. The x-component of 
momentum that enters the x-face of the element at point x through 
20 convection is Az*Ay«{p»v x «v x } x , and the rate at which momentum 

leaves the element at point x+Ax is Az*Ay*{p*v x «v x } x+Ax . The rate 
at which momentum enters the element through the y-face is 
Az*Ax*{p*v»v__} v , and so on. Adding all six components of momentum 
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in the x direction yields, 

Ay Az [ (p v x v x )\ x - (p v x v x )\ x ^ ] ♦ 

7.6 

AxAz [ (p v y v x )\ y - (p v y vjl y+Ay ] + 
Ax Ay [ (p v z v x )\ z - (p v z v x )l^ A , ] 

The viscous stress induced x-component of the normal 
deviatoric stress at point x is Ay«Az*{T' xx } x and likewise for the 
other stress components. Summing all these stresses acting in the 
x direction produces, 

Ay Az IT*„\ X - r^l^, ] + 

7.7 

Ax Az [T*„\ y ~ r^JyAy 1 + 
Ax Ay [ T*„l a - T f zx \ z . Lz ] 

In addition to the deviatoric stresses there is also the net 
hydrostatic pressure, P*, acting on this volume of material, and 
the forces exerted by gravity, g. Summing the x-component of 
these forces yields, 



Ay Az [ P* l x - P* l x+Ax ] + Ax Az A y [ p g x ] 

7.8 



Lastly, the total rate of x-momentum accumulation within the 
element is simply, 



(Note that: Ax/At = v_) 



218 



Inserting these equations into the momentum balance 
equation, dividing by Ax»Ay»Az, and taking the limits as At, Ax, 
Ay, and Az go to zero produces the x-component of the equation of 
motion, 

d 9V X _ _ ( d pv x v x + apVyVx + d pv z v x 
a t ~ d x d y d z 7. 10 

, a a* a a r'^ a p« 

" { inf + "aT* + "TIT 17 + p9r * 

Similar results will be obtained for the z- and y-components of 
momentum and force yielding, 



3t Jx d y d z 



7.11 



3x 3y 3z 3y ry 



IP 
IB 

- and 



3 py, _ ( 3 pv z v x + dpv y v z + a py^ } 
10 i» 3 ' 8X 3y 3Z 7.12 



I* 1 6x 3y 3z ' az 



These three equations may be rewritten as, 



iL£Z . -V(ovv) - VP* - V r' ♦ pG 
at 7.13 



where V is the velocity vector, 6 is the gravity vector, and T' 
is the deviatoric stress tensor. The gravity vector G has 3 
15 component terms. The stress tensor T # has nine component terms, 
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as does the momentum flux term, pW. 

Equation 7.10 may also be rearranged into, 

d PVx + / 9 f> V x V * + d P v y V x + d P v s v x ) _ 
a t d x d y d z 

1 3x 3y d z ' d x py * 

where the left hand side of this equation is, 



7.14 



3t 1 3* d y d z 7.15 



i3 These four terms may in turn be expanded into the expression, 
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3 v r d v x dv x d v x 

K a t r * d x r y d y r z d z 

m + pvv —5 — — + pv x -r-X- + pV x - /.XO 

|3 .* K dt*dx y dy ' d z ' 

SI 

S ! J Refer back to equation 7,2. Substituting this expression 

13 into the last term of equation 7.16 yields, 

|£ 

d v x d v x dv x d v x 

+ P v x + P v y -5-7; + P v « "~ ~" 



10 y d t v * x d x dy rw * d z 7.17 

Recall that the definition of the 'substantiative' or material 
derivative, Dc/Dt, is, 

DC _ dC + dC + dC + ac 



Dt at * ax y ay *dz 7. is 
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# • 

where C is a dummy variable. From this definition equation 7.17 
may be rewritten as, 



D v, d V- d v x dv x d v x 

-DT - P ^f + Pv *-df + pV * 37 + PV '^T 7.19 



Expressions similar to equation 7.19 may be derived for equations 
7.11 and 7.12 and equation 7.13 may then be rewritten as, 

K d t r 7.20 



Note that the two stress terms, the net hydrostatic pressure P*, 
and the deviatoric stresses associated with the viscous forces, 
T', are expressed separately. 

7.3 The Governing Constitutive Equations 

To proceed any further in the solution of equation 7.20 it 
is necessary to define the appropriate constitutive response for 
both the viscous stress term and the hydrostatic stress term. 
Let's restrict interest to the simple Newtonian fluids. The 
Navier-Poisson constitutive equation for Newtonian fluids was 
derived on the assumption that the viscous forces acted not only 
deviatorically, but also hydrostatically . This assumption 
produces a set of stress equations of the form, 
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which are analogous in form to the equations of isotropic linear 
elastic solids. The total stress T 1 ^- is defined to be the sum of 
the applied stresses, T ij# and the thermodynamic stress p. 
Equation 7.21 may be expressed equivalent ly as, 



T ij + p 8 ij - c ijki D ki 



7.22 



or 



Zjj - ' P &ij + Cijkl D kl 



7.23 



?0 

The normal stresses, the T n , may always be written in terms 



jn of a deviatoric component, T', and a hydrostatic component, T h , 
10 la i.e., T i4 = T' n + T h iS . The hydrostatic stress component is also 
equivalent to the total mean normal stress or pressure, -P, 



P 



rryh Tfcfc ^11 + ^22 + ^33 _ 

* u " " ^ 



IU r « " 1 3 • * 7.24 

in 
p 

where the (-) sign indicates compression. From this definition 
and + p = -P + p, the net applied hydrostatic pressure is 

15 then simply -a^ = -P* = -P + p. Obviously, if a lM = P* = 0, then 
P = p. Equation 7.24 leads to the generalized stress expression, 



Demanding isotropy (and following the development of a similar 
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set of equations for the linear elastic solids) equation 7.23 
then becomes, 



13 



7.26 



where X and M are the Lame' coefficients, analogous to similar 
coefficients in the theory of elasticity. This is the Navier- 
Poisson equation for isotropic Newtonian fluid flow. The C lJki 
coefficients of equation 7.23 are now defined as, 



X + 2\i 
X 
X 

0 
0 
0 



X 

X + 2\i 
X 

0 
0 
0 



X 
X 

X + 2ji 

0 
0 
0 



0 
0 
0 

0 
0 



0 
0 
0 
0 

0 



0 
0 
0 
0 
0 

2*1 



7.27 



where M is the dynamic viscosity. This matrix is identical in 
form to that for linear isotropic elastic solids. 

The normal rate of deformation terms, the D^-, may also be 
written in terms of a deviatoric rate of deformation, D' n , and < 
mean volumetric rate of deformation, D klc /3, i.e., 



v • v 



- D f 



ii 



D X1 + D 22 + £> 33 



7.28 



ii 



This function in turn leads to a generalized expression for the 
rate of deformation, 




where D|j is defined by equation 5.3. 

Substituting equations 7.25 and 7.29 into equation 7.26 to 
separate deviatoric and hydrostatic responses yields , 

J 7.30 



P This equation may now be separated into an equation that defines 
=0 the hydrostatic response, 



7.31 



iff , de 

y P - p - (X * iiiJOj* - P - JCD tt - P - K-^f 

IB 

W -> X" - A, + -|n 

□ 

N 

*~ and another that defines the deviatoric response, 
10 1J P iJ 7.32 



Equation 7.31 now has two additive constitutive functions 
that define the total constitutive response under hydrostatic 
loading, the thermodynamic equation of state p, and a viscous 
term, K*de v /dt. e v is the volumetric deformation and K is the 
15 'bulk viscosity'. Equation 7.31 results from the assumption of a 
viscous constitutive response for hydrostatic loading identical 
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to that for deviatoric loading. This equation is incorrect, for 
it implies that some net sudden (shock) impulse hydrostatic 
stress, defined as a i{ = P* = P - p at time t = 0+, will produce 
only a viscous response, i.e., 

p*-p-p--u + 4^ D k* - - KD kk - - * — 



There is simply no real material that behaves in such a 
fashion. Consider a volume of gas that is subjected to a sudden 
shock impulse stress generated by a piston, an explosion, or 
P whatever. Equation 7.33 implies that the gas molecules, which are 
10 disconnected and floating in free space, do not recoil 

iff elastically from the impact of other gas molecules, or the impact 

W 

ig of a solid interface, the surface of our piston, etc. Hence, 

ig 

. there is no shock wave or any physical mechanism to propagate a 

H 

;% shock wave, or any elastic (acoustic) stress wave for that 
15 matter. Furthermore, at the atomic level equation 7.33 also 

}~ implies that the electrostatic repulse forces operating between 
(at the very least) the nuclei of colliding gas molecules, or the 
gas molecules and the surface of the piston, have no influence on 
material behavior. This result is clearly physical nonsense and 
20 comes about from the assumed form of the governing constitutive 

function, equation 7.21. This result is also why the equations of 
state, and not the Navier-Poisson equation, are used for wave 
propagation studies. 

Note also that equation 7.31 is analogous in form the 
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Maxwellian viscoelastic fluid model, where this hypothetical 
constitutive function is the sum of an elastic response term (the 
thermodynamic equation of state, p) and a viscous response term 
(K*de v /at) . No known material, fluid or otherwise, behaves in 
this fashion for the viscous term allows for unconstrained 
deformation for any level of applied stress. 

So, the correct form of the complete constitutive equation 
of Newtonian fluids must take the form of two independent and 
linearly additive responses and may be written as, 

where the volumetric pressure term, call it P = P(P*/P#0#p/e v , t) , 
and the deviatoric term are characterized by two different and 
uncoupled functions. P = P(P*,p,0/P, e v ,t) is an equation of state 
that defines volumetric deformation, e v , as a function of net 
stress hydrostatic stress P* = P - p, temperature 0, density p, 
and time t. 

Note that equation 7.34 is analogous to the 'Navier-Poisson 
equation of Newtonian fluids with no bulk viscosity'. The form of 
this function is not new (except for the recognition of 
viscoelastic behavior for the equation of state and the required 
uncoupling between the two responses) . Historically this function 
has been derived by imposing the approximation that bulk 
viscosity may be ignored for most applications. This is 
equivalent to demanding X » -2/3/x in equation 7.30, the so-called 



Stokes condition. In our case, equation 7.34 represents the 
complete and proper solution, not an approximation. 

7.4 Final Solution for the Equations of Compressible Plow 

Now that the proper form of the constitutive equation has 
been determined, we need only substitute the appropriate 
expressions into equations 7.20 to obtain the final form of the 
equation of motion for compressible fluids. Equations 7.3, 7.25 
and 7.32 yield the following solution for deviatoric stress, 

m/ 



10 *U J 7.35 

i '"■ 

in 

y 

10 From equations 7.26 and 7.28 the expanded forms of T' fj and T' M 



10 



P 



are then, 



SJ T'u - -2|i ■& + 4l»'< V * v > 



|t 5 

in 

13 

Isfe 



** 3 7.36 



Symmetry in the shear stresses is guaranteed by symmetry in D^- 
15 (note the sign convention change to conform to the sign 

convention of equation 7.20, which comes about from equations 7.6 
- 7.8) . 

Now, expanding the x-component of equation 7.20 yields, 
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Dv 



D t 



? - - ( 



d x 



dy 



a z 



) - 



d p* 

d x 



?9 X 



7.37 



Substituting the appropriate expression for stress, as defined by 
equation 7.36, produces, 



7.38 



Similar expressions will be obtained for the y-component, 

- -|t t p ?y ♦ » - 4 p( V-v>] + 

Dt dy r y dy 6 y 3 r 7.39 



in -k^t-df + -33T )] + "3* 14 ~5"y "3lf )] 

hi 

ifl 

A " and the z -component, 

Q 



7.40 



13 a [M 4^ + + a^^*^,] 

|4 dx r d x dz dy r d z dy 

These three equations, coupled with the temperature, density, and 
10 a viscoelastic equation of state (bulk compressibility function) , 
P = P(P*/P/0#P/e y# t) , the temperature and density dependent 
viscosity function /i = M(0#/>), and the initial boundary 
conditions, completely determine the pressure, density, and 
velocity components of a flowing compressible isotropic Newtonian 
15 fluid. Assuming a constant viscosity i.e., invariant 0 and p, 
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these three equations may be rewritten in a more compact format 
as, 

V Dt a. 3 7.41 



where the net hydrostatic pressure function P* is defined in 
5 terms of the equation of state, P = P(P*#p,e v , t) , which is a 
viscoelastic function. 

Equations 7.38 - 7.41 may also be obtained by substituting 
equation 7.34 directly into Cauchy's equations of motion, i.e., 

fi „ DV 3 DV 1 

g v-r ♦ pa - p% -> ♦ P9J - p-^ . 7 

(n 

10 |Q Note that equation 7.20 is equivalent to equation 7.42, only 

in 

equation 7.20 has a reverse sign convention for stress. 

13 

\z 7.5 Prior Derivations and the stokes condition 

ill 

*** The Complete Navier-Stokes Equations 

Assuming for the moment that there is such a property as 
15 'bulk viscosity' , equations 7.20 and 7.26 may be used to derived 
the complete traditional form of the Navier-Stokes equations for 
compressible fluids. While the viscous deviatoric stress terms 
are unchanged by this assumption, i.e., 
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The hydrostatic viscous stress terms now become, 

" p 3x 4 7.44 

Substituting these equations into equation 7.37 yields the 
x-component solution, 

p Similar expressions will be obtained for the y-component, 



u t dy y dy 3y 



7.46 



dx r d x dy dz r d y dz 

m 

13 

and the z-component, 



IU 



7.47 



D t dz z dz d z 

3x 3^ dz dy r d z dy 



10 Letting X and M be constant, these three equations may be 
rewritten as, 



y Dt r 7.48 



Equations 7.45 - 7.48 are the complete form of the Navier-stokes 
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equations for compressible fluid flow. However, these equations 
are based upon an assumed form for a governing constitutive 
function, equation 7.21, that in part describes unreal material 
behavior, so they cannot be correct. 

The Stokes Condition 

Return to equation 7.30 and impose the Stokes condition, 
i.e., the assumption that viscous stress effects are negligible 
for volumetric deformation. This leads to the relationship 
X = -2/3/4/ equivalent to K ~ 0. The only viscous effects allowed 
are those associated with deviatoric loading and deformation. The 
Stokes condition now eliminates the assumption in equation 7.21 - 
7.23 that hydrostatic loading and the associated volumetric rate 
of deformation are also governed by the same basic constitutive 
relationship of Newtonian fluid flow that governs the deviatoric 
rate of deformation. 

Substituting the Stokes condition term X = -2/3/i into 
equation 7.30 produces K = X + 2/3/i = 0 and, 

a result similar to equation 7.34. This equation likewise implies 
20 that there are two different, and uncoupled, constitutive 

functions now governing fluid deformation: the thermodynamic 
equation of state, p, and the constitutive equation relating the 
deviatoric stresses to the deviatoric rate of def ormation, 
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2/* # D'jj. Substituting equation 7.49 into the governing equations 
of motion 7.20, letting p = P, the result will again be equations 
7.38 - 7.41. 



O 

in 

l>i 

i.SJ 
!» 

IS 
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SECTION 8 

APPLICATION EXAMPLE - NUMERICAL SOLUTIONS OF THE LINEAR 
HEREDITARY INTEGRALS 

The hereditary integrals provide a means of obtaining 
5 closed-form solutions for the constitutive equations under 
dynamic loading or deformation conditions. For steady state, 
periodic loading these integrals degenerate into a unique format. 
However , not all loading is periodic, nor is it always steady 
state. Furthermore, not all applications (e.g., beam, plate, 
10 q shell, etc. theories) lend themselves to closed-form mathematical 
.p solutions for all possible boundary and loading conditions. In 
m addition, static loading of a body of material may produce 

: r : 

-!g hereditary effects because internal stresses can vary as a 
^ function of deformation and time. 
15 P As a consequence of constraints inherent to any closed-form 

;m 

lU solution, including the various application theories, the 

in ■ 

St * 

13 hereditary integral functions need to be adapted for numerical 
applications, such as finite difference or finite element 
methods. These numerical methods, which incorporate material 

20 constitutive properties, can provide accurate (albeit 

approximate) solutions to a wide variety of engineering and 
scientific problems. A scheme is therefore needed for solving the 
various forms of the hereditary integrals numerically, for both 
static and dynamic loading or deformation of all forms. We also 

25 need to be able to invert the hereditary integrals. Once 
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developed this method may then be incorporated into the higher 
order finite element , finite difference, etc. numerical solution 
schemes, thus forming the technical foundation for all derivative 
forms of computer software. 



8.1 The Linear Constitutive Equations 

Let's restrict interest to the strain limited linear 
approximations of viscoelastic behavior where constitutive 
properties are given by the functions, 



•? 
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10 in and 
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IB o £ (t) - [a ij ]f E e j + [b ij ]f D (t)e j 
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8.2 



which define response to a static load or static strain. Recall 
that these equations may be rewritten as, 



• lit) - t k iJ (t)]f*(t)o j o^t) - [ l iJ (t)]f(t)e J 

8*3 



15 



or 



«i(t) - [D ij {t)]o j 



o^t) - [ C„(fc)]e 



8.4 



where 
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+ n u f* D (t) 



f* s + f* D (t) 



f+{ t ) - £+* + /♦^(t) 



8.5 



and 



I } f* B + [ ^ ] f* D (t) 



ii 



8.6 



and similarly for IjjCt), f(t), and C-jCt). Note that the k^-tt), 
Dij(t), etc. are not symmetric in time unless each static a } and 
€j is applied simultaneously. They are however, always symmetric 
P at time t = 0+ and t = «. Equations 8.3 may readily be inverted 
;Q to yield solutions for the static stresses, the o- 9 or the static 
jfj strains, the €j, i.e., 



to 

\U oj - e^t) [ kyit)]- 1 f*(t)- x ej - a^t) [ l^it)]' 1 f(t)' 1 

1 ^ c 

i, 5 

j~j From previous discussion we know that the hereditary 

^ integrals of equations 8.1 and 8.2, which quantify stress and 
strain under dynamic loading or deformation, may be expressed as, 



e^t) - ^f.'^lt) + n ij jf* D (t-t')do J (t') 



8.8 



15 



and 



Oi(t) 



'1J 



f^jlt) + b u f fit - t')<fej(£') 



8.9 
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or similarly as, 

e^t) - m u f**o.,(t) ♦ n u J f* D (t - t') ^j}?* dt 



dt' 8 . 10 



and 



r^t) - a^fejit) * b tJ [ f D (t - t')^p-dt' 



dt' 8.11 



5 or, after some manipulation, 

□ «< (t) ■ * **> l[y (t,) ^dlt'.'t^ dt ' 8.12 



and 



{ (t) - aijftejU) ♦ b„ / €j(tO a^/fc/ 



IU 



dit - t') 8.13 



Recall that equations 8.1 -> 8.6 are limiting solutions to the 



10 :^ hereditary integrals for static loading or deformation. 

8.2 Incrementing the Constitutive Functions 

Figure 8.1a illustrates some time dependent stress function 
a(t) which produces the time dependent delayed elasticity 
response f*°(t) (Figure 8.1b). f* D (t) may be approximated as the 
15 sum of a series of incremental responses (Figure 8.2a) resulting 
from the application a series of incrementally applied stresses, 
the 6 n a(t n ) (Figure 8.3a). The better the stress function 
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approximation (Figure 8.3b) the more accurate the strain history 
approximation (Figure 8.2b). Differentially small stress 
intervals produce the hereditary integral functions of equations 
8.8 and 8.9, as well as other various equivalent forms of those 
functions, e.g., equations 8.10 through 8.13. 

Let's review the method of establishing the hereditary 
integral function 8.8, but refrain from converting the resulting 
algebraic equation into an integral function. Consider equation 
8.1 which may be rewritten as, 



If! where, for example, e^t) is simply, 

ui 
m 

ifl eil (t) - €* n + e D (t) 11 - ra^ f* B o x + n 11 f* D (t)a 1 



8.15 



'-4 



in 

13 



By choice of definition m n = n n = 1. Figure 8.4 is the graphical 
illustration of the two strain components. 
15 Let's now incrementally increase the stress with time. 

The time interval between incremental stress increases being 
irrelevant. Each incremental increase in stress, the S n a^ 9 
produces linearly additive impulse and linearly additive delayed 
elasticity strain responses (Figure 8.5). The solution for the 
20 elastic impulse response as a function of time is then, 
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«*ii - ^lE^^n^It.) - '**E*»°i< t J 

"-° n- ° 8.16 

- ^ f+* 6 0 o 1 (0+) + ai ll f* B 6 1 a 1 {t 1 ) + ... 
. . . + m ll f* B & D o 1 (t a ) 

where each S n a^(t n ) represents a constant valued incremental 
increase in a t applied at some time t n . t 0 = 0+ by choice of 
definition. The corresponding delayed elasticity response will 
likewise be the sum of the incremental responses, 

a-0 8 . 17 

- n lx f* D (t) 6^,(0 + ) + n lx f* D (t-t 1 )b x o 1 {t 1 ) + 
... + n xl f* D (t-t n ) b n a 1 (t a ) 
where the solution for e n (t) = € E U + c^^t) is, 

Similarly, the solution for the six € f (t) now becomes, 

«i<« - ai^E'** 8 + ^E^'^-g 8 ^'^ 8#19 

where the only applied stress is a v 

For convenience the full three dimensional form of equation 
8.19 may be simplified by requiring that incremental increases in 
the 0j occur at the same interval in time, t n , although the 
magnitude of the incremental increases, and the number of 
increases in stress may be different for each o-. Equation 8.19 
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may then be expanded to, 

The time interval t n is incremented from the time of application 
of the first incremental stress component to the last increment, 
5 regardless of which c- } is applied first and which is applied 

last. There are nj increases in stress for each Oy Each stress 
may also have a different t 0 as well as a different time of 
application of the last incremental increase in stress. 
13 It is from equation 8.20 that equation 8.8 is obtained. The 

10 \Q conversion is made by simply driving each finite incremental 
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stress to a differentially small value and converting the 
summation to an integral, which is simply the sum of all the 
differential component strains through the time interval t. The 
5^ integral is evaluated by integrating over each increment of 
15 }~ stress using the proper t n in the argument of f * D (t-t n ) . Equation 
8.10 is also obtained directly from equation 8.20, where the 
substitution da(t') = (da(t') /dt') # dt' has been made for da(t'). 



in 
p 



Alternative Solutions for The Creep Compliance Functions 

Each strain component is linearly additive, despite any non- 
20 linearity in f* D (t), so we may rewrite equation 8.18 as, 
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[in^ f* s i^it^ + /3 1X f* i> (t-t 1 )6 1 o 1 (t 1 )] + 
or more compactly, 

where the total strain at time t is simply the sum of the 
incremental strains. Equation 8.19 then becomes, 

which then leads to the following equivalent of equation 8.20, 



8.21 



8.22 



Let's now rewrite equation 8.21 in a form similar to that of 
equations 8.3 and 8.5, i.e., 



, . r jn ll f* B ♦ i3 lx t) 



which is equivalent to, 
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en(t) - E *ii<t-t n ) f*(t-C„) J^o^t,,) g>26 



where 

m„ f* B + iJ„ f* fl (t-tj 



/** + £* J '(t-t n ) 8.27 



and f*(t-t n ) = f* E + f*°(t-t n ). So, the value of k n (t-t n ) varies 
5 in time for each incremental value in stress. For example, when 
time 0 > t > t, then k^Ct-t,,) = k^tt) because t 0 = 0+, and is 
q equivalent to, 



/** + f* D (t) 8.28 

Ml 



w 

03 



When time t t > t > t 2 then k^t-t,,) = k^tt-t,), i.e., 



N ..... ^ f«*-m u f^tt-t,) 

10 IU k^t-tJ f+t*Ht-t,) 8.29 

in 
a 
\* 

and so on for all t n . Equation 8.23 may now be rewritten as, 

Mt) - 5>i<t„) " E *n(t-t a ) f*(t-t„) fi.uja^t.) 



8.30 



where the k n (t-t n ) are given by, 

f* B + f* D (t-t n ) 



^i(t-t a ) 



/♦* + f* B (t-t n ) 8.31 
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Given that each incremental increase in the occurs at the same 
point in time t n , the delayed elasticity constitutive function 
f*°(t-t n ) has the same value for each incremental $ n(1) <7j* So, the 
three dimensional form of equation 8.30 may be rewritten as, 



- £ kuU-t n ) f*(t-t n ) b nU) o 3 (t nij) ) 



8.32 



where the k^Ct-tJ are, 

ju ( t t ) m ** f * S+n U f* D i*-*J 

K V [C - C " } f* + £*Ht-t a ) 8.33 

%B Note that when the approximation m^- = is applicable then 

jfj k|j(t-t n ) » m^-. Also note that equation 8.32 is identical in form 

ill 

10 ijg to Hooke's law at any time t. By letting each incremental 

iB 

increase in stress, * n(J -/7j# be applied at the same incremental 
time t n , we also insure that the incremental equation 8.33 



iU 

in 

13 



produces a symmetric matrix of time dependent coefficients, the 
incremental k f j(t-t n ). 



15 The Relaxation Functions 

Following the above discussions equivalent solutions for the 
relaxation response functions are given by, 



where the ljj(t-t n ) are defined as, 
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l)K " } f*+f°(t-t D ) 8.35 



and f(t-tj = f E + f D (t-t„). These functions are likewise 
identical in form to Hooke's law at any time t. 

8.3 Inversion of the Hereditary Response 

From the above solutions for creep compliance we may obtain 
solutions for the loading induced stress histories if the six 
€ i (t) are known as a function of time. As we have seen, all 
P solutions for strain are linearly additive. Each component of 

■q 

kU stress is also linearly additive where the value for each of the 



10 in tfj is simply, 

IS 



°-o 8.36 
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Let's now presume that some body of material has been 
H subjected to loading functions, producing about a small finite 

; • 

element of that body the surface tractions Cj(t) . Let's further 
15 presume that, aside from the material properties, i.e., the m-. 

and n ij# f* E , f D (t), etc., that we know the strain histories €.(t), 
which may obtained from the displacement histories. The variation 
of the Ojit) as a function of time (and therefore also strain) 
are unknown and it is these stress histories that we seek to 
20 determine. An example of such a case would be a structure 

subjected to one or more loading functions that produce bending, 
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twisting, displacement, etc., where the € f (t) are determined from 
the displacement histories, strain gauges, etc. 

As discussed, the time dependent and presumably (albeit not 
necessarily) continuous stress functions may be approximated as 
step functions which increase in time in very small but finite 
increments. Let the initial applied incremental stresses be the 
*o<j> a j* 0+ *' wh i ch are a11 unknown. At any time t > 0+ the e^t) 
are known. Let t t be some small time interval greater than 
t = 0+, the instant the first incremental stress was applied. 
From equations 8.32 and 8.33 the solutions for the associated 
strains at time t 1 are then simply, 



where 



my £* B + n ±j f* D {t x ) 

f** + f* D {t x ) 8.38 



Given that the tf^Ct^ are known from physical measurements the 
solutions for the incremental stresses, the * 0 <j) a j(° + )' at time t t 
are then, 



where these were applied at time t G = t = 0+. 

At time t,+ let another constant increment in the stresses 



be applied, the f 1(J) a j (t 1 'f) , which are unknown. At time t 2 > t t + 
the solution for strain from equation 8.32 is, 



e^tj) - t k^it,) ] f*(t 2 ) « o(J) o.,(0-O * 
where at time t 2 , 



8.40 



*« (t » } f** + f**(t: a > 8.41 



and 

§3 

J t.* * tl'itt-tj 8.42 

CSS 

in 

*f 5 

Id 

Equation 8.40 may also be rewritten as, 



iU 

in 



€i(t 2 ) - e i(1) (t 2 ) + e^j^ta-ti) 



8.43 



10 The values of the initial constant valued S QiJ} c } lQ+) were 



determined from equation 8.39. Equations 8.36 and 8.37 then yield 
solutions for the * i(0) (t) which continue to increase in value 
with increasing time unless the asymptotic limit is reached for 
f* D (t). Knowing the values of the e^t) at times t, and t 2 the 
15 solutions for the incremental * i(2) (t 2 ) are simply , 

•iwt^-^i) - € i<t 2 ) " «iu>< t 2> q AA 

8 • 44 
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From equation 8.40 the solutions for the incremental 'icj> a j Ct t ) 
are then, 

*iu>°><*i> " le^t 2 ) - e i{1) (t 2 ))[k iJ (t 2 -t 1 )]- 1 f+itz-tj- 1 8.45 
So, for time t 2 the approximate total stress solutions are, 



which holds true until the incremental time t 2 + when the next 
Q unknown incremental stresses were applied. 
<p Let's iterate once more, after which we may write the 



\f\ generalized solutions. At time t 2 + the unknown incremental 

* - 3 

10 j| stresses * 2 <j> a j 1*2* are applied. The next application of 



incremental stress was at time t,+ , so we want to evaluate for 

a 

P. strain at time t,. From equation 8.32 the solution for strain in 
terms of the incremental stresses is, 



e ± (t 3 ) - I k 1:l (t 3 ) ] f*(t 3 ) 8 o(i) a i (0+) + 



8.47 



15 where at time t 3 , 



. . , f* s + f» p (t,-t 2 ) 



kjjttj) and kjjCtj-tj) are given by equations 8.41 and 8.42 
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evaluated at time instead of time t 2 . Equation 8.47 may now be 
rewritten as, 



. e^tj) - e iU) (t 3 ) + e i(2) (t 3 ) + e i(3) (t 3 ) 

8.49 



The constant valued * O (j> a j (0+ ) and S MD a i^ t 0 were determined 
by equations 8.39 and 8.45. The contributions of the € i(1) (t) and 
the € i(2) (t) continue to increase in value with increasing time. 
Knowing the ^(tj) the solutions for the € 1(3) (t 3 ) are, 



e i(3) (t 3 ) - e^tj) - € i(1) (t 3 ) - c i(2) (t 3 ) 
13 8.50 

*□ 

M 

ifi The corresponding solutions for the 5 2(j) aj(t 2 ) are then, 

id 



ft 2(j »^(t 2 ) - € i(3) (t 3 ) tJc i:/ (t3-t 2 )]- 1 £*(t 3 -t a ) 



-1 



10 m 8.51 

- te^t,) -e i(1) (t 3 ) - c i(2) (t,)] [* ii (t 3 -t a )]- 1 /♦(tj-ta)- 1 

Hi 

^ Lastly, for time t 9 > t > t 3 the total stress solutions are, 

in 

IS 

- a jit) - lowOjiO*) + ft^o^) + 8 2(i) a i (t a ) 



From the above examples it should be clear that for the 
finite time interval the solutions for the incremental 



15 stresses applied at time t n are, 
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* a u>« j < t n ) - e, (ml) ( t) [Jc^ ( fc ml - t n ) ] - 1 /• ( t^- t a ) - 1 8 . 53 

- [« 4 (t^i) " E'WWl [^(t^-tj]" 1 ^(t^-tj- 1 

n-0 

a solution still identical in form to Hooke's law. From equation 
8.52 the total stress solutions are given by equation 8.36. 

Similar functions may readily be derived for the relaxation 
5 functions of equation 8.34 and 8.35. Those derivations are now 

straight-forward and warrant no further discussion. Obviously we 
need not confine interest to numerical solution of the hereditary 
u integral functions of the form given by equations 8.8 -> 8.11. 
=p Equations 8.12 and 8.13 can be reformulated in a similar fashion. 
10 m In the end however, the same basic procedures must be followed 

E . 3 

]q for obtaining solutions for the hereditary responses, and well as 

^ the inversion of those responses. 

13 

— 8.4 Additional Comments 

P As mentioned previously, we have chosen to increment each 

15 stress function at identical t n intervals. The reason for this 

choice being that the resulting k^ft-t^ and l--(t-t n ) are always 
symmetric. Because of this symmetry the inverse of these 
matrices, the kjjit-t^" 1 and l^-tt-t^" 1 are also symmetric. 
Furthermore, the constitutive equations may be written in a form 
20 identical to Hooke's law so we automatically know the solutions 
for the inverse relationships of k S j(t-t n ) and l-j(t-t n ), 
regardless of the degree of material symmetry. For example, for 
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isotropic materials Hooke's law may be written in the form of, 



- [ D u } f**o 



8.54 



and the inversion as, 



8.55 



where f E = l/f* E . For isotropic materials f E = E, which is Young's 
modulus. The solution for the coefficient matrix D Vl is, 



1 


^12 


^12 


0 


0 


0 




1 




0 


0 
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0 


0 


0 
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0 












0 
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The material coefficients D 1V D 12 , and D u are defined as, 



A, - 



a + P 



P(3a + 2P) 



- 1 



8.57 



*12 



2p(3a + 20) 



. - v 



^44- 2(^-^2)- -j " 2(1 + V) 



where i; is Poisson's ratio, a and /? are the Lame' coefficients 
defined as, 



a - 



(1 + v) (1 - 2v) 



P 



2 (1 + v) 



8.58 



From Hooke's law [D^]" 1 is, 
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a+2P a 
a+2f} 



where 



a 
a 

a+2p 



0 
0 
0 

p 



0 
0 
0 
0 

p 



8.59 



a + 2P - 



1 - v 



(1 + v) (1 - 2v) 



8.60 



So, by simple extension the solution for the k^t) for 
isotropic materials is, 



hi 

i 



!U 
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Sea 



k 12 (t) 
1 



* 12 (t) 

*1 2 U) 
1 



0 
0 
0 

2(l-k 12 (t)) 



0 0 
0 0 
0 0 
0 0 
2(1-Jc 12 (t)) 0 

2(l-ir 12 (t)) 



8.61 



and the solution for l/k^-ft) is given by equation 8.59 where, 

-k 12 (t) 



[l - ic 12 (t)] [l + 2Jc 12 (t)] 



P 



2(1 ♦ k 12 (t)) 



8.62 



a + 2P - 



1 + * 12 (t) 



[1 - Jc X2 (t)] [1 + 2* 12 (t)] 



Thermal Effects 

10 In general, the equations for thermodynamic expansion and 

contraction may be developed in the same manner as the 
constitutive functions. However, there is one important caveat. 
Temperature effects operate like a hydrostatic load and there are 
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no shear equivalents. So, if a material has directional 
singularities to it's constitutive properties, these 
singularities will dominate the thermal properties in the axial 
direction of the singularity. 

8.5 Discussion of Applications 

The simplest possible form of the constitutive equations for 
linear viscoelastic materials has been unknown. In general, the 
basic fundamental formulation for the constitutive functions has 
been assumed to be of the form, 



n 



10 *□ 8.63 



m where the D s; (t) and C iS (t) have been characterized as separate 



and independent functions, in a fashion analogous to linear 
elasticity. As we have seen, this assumed independence is a 
violation of the fundamental physics of material deformation and 



SU 

in 

15 H the fundamentals of strain theory. From earlier discussions we 



know that the D^it) and C rj (t) have solutions of the form, 

* 3 8.64 

C^it) - [ a ±J ]f B + t b £j ]f D (t) 

and are not symmetric in time unless all loads have been applied 
simultaneously . 

20 Numerous numerical methods and associated computer software 

have been developed based upon equations 8.63, and associated 
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assumptions regarding material behavior. These methods may be 
readily adapted and modified for the proper form of the 
constitutive eguations and material properties through the use 
equations 8.64 and the methods for solving and inverting the 
hereditary integrals. 
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